Sharedcloud-examples / sage / 2018-02-15-162612_Picard_Numérico_Integral_y_ODE.sagewsOpen in CoCalc
Iteraciones de Picard Numéricas por Cuadratura y por solución numérica de ODE
# -*- coding: utf-8 -*-
# Iteraciones de Picard Numéricas por Cuadratura y por solución numérica de ODE.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from sage.plot.colors import rainbow
# Picard iteration using quad (Cuadratura)
def picard1(f,t,y,N):
    TOL = 1e-8; # default, reset if desired
    z = y[0]*np.ones(t.size)
    for k in range(1,N):
        I=quad(fint,t[k-1],t[k],args=(f,t,y),epsabs=TOL,epsrel=TOL)
        z[k] = z[k-1]+I[0]
    y = z
    return t,y
#-------------------------------------------------------
# Picard iteration using odeint (Solución numérica de ODE)
def picard2(f,t,y):
    TOL = 1e-8; # default, reset if desired
    y=odeint(fint1,y[0],t,args=(f,t,y),rtol=TOL,atol=TOL)
    return t,y
#-------------------------------------------------------
# ODE rhs interpolated at t from (told,yold)
def fint(t,f,told,yold):
    spline=interp1d(told,yold,kind='cubic')
    dydtint = f(t,spline(t))
    return dydtint
def fint1(y,t,f,told,yold):
    spline=interp1d(told,yold,kind='cubic',fill_value='extrapolate')
    dydtint = f(t,spline(t))
    return dydtint
#---------------------------------------------------
def f(t,y): # ODE
    dydt = np.exp(y+t)
    return dydt
def f1(y,t):
    dydt=f(t,y)
    return dydt
#---------------------------------------------------
# initialization for Picard iteration.
def init(t,y0):
    y = y0*np.ones(t.size) # standard. If modified, make sure the new function satisfies the IC
    return y
#--------------------------------------------------------
# plots the direction field for the scalar ODE y'=f(t,y)
# in domain specified by window using a 21x21 grid of arrows
# use: dirfield(@f,[tmin,tmax,ymin,ymax])

def dirfield(f,window):
    t = np.linspace(window[0],window[1],21)
    y = np.linspace(window[2],window[3],21)
    T,Y = np.meshgrid(t,y)
    k=window[3]/window[1]
    # evaluate slope and normalize
    S = f(T,Y);
    #N = np.sqrt(1+S**S); dT = 1/N; dY = S/N      # Se omiten estas instrucciones de normalización ya que quiver se encarga
    #del reescalamiento de los vectores.
    Q=plt.quiver(T,Y,1,S/k,units='width',scale=100)
    grd = plt.grid(True)
    plt.tight_layout()

#-------------------------------------------------------------------
# Iteraciones de Picard Numéricas por Cuadratura
tspan = [-0.5,0.296]; yspan = [0,5]; y0 = .3; N = 10;c=rainbow(N+1)
window = [tspan[0],tspan[1],yspan[0],yspan[1]]; #t = tspan
npts = 100; # #pts for integration
t = np.linspace(tspan[0],tspan[1],npts)
y = init(t,y0)
p=plt.plot(t,y,'r-',linewidth=2.0)
for n in range(1,N):
    (t,y) = picard1(f,t,y,npts)
    p+=plt.plot(t,y,color=c[n],linewidth=2.0)
dirfield(f,window); # direction field, see Lab1
ye=odeint(f1,y0,t); # "exact" solution
p+=plt.plot(t,ye,'-',color='gray',linewidth=2.0); # plot "exact" solution
plt.ylim(yspan)
plt.show()



(0, 5)
#Iteraciones de Picard por solución numérica de ODE
tspan = [-0.5,0.296]; yspan = [0,5]; y0 = .3; N = 10;c=rainbow(N+1)
window = [tspan[0],tspan[1],yspan[0],yspan[1]]; #t = tspan
t = np.linspace(tspan[0],tspan[1],100)
teval = np.linspace(tspan[0],tspan[1],200)
y = init(t,y0)
p=plt.plot(t,y,'r-',linewidth=2.0)
for n in range(1,N):
    (t,y1) = picard2(f,t,y)
    y=y1[:,0]
    spline=interp1d(t,y,kind='cubic')
    yint=spline(teval)
    p+=plt.plot(teval,yint,color=c[n],linewidth=2.0)
dirfield(f,window); # direction field, see Lab1
ye=odeint(f1,y0,t); # "exact" solution
p+=plt.plot(t,ye,'-',color='gray',linewidth=2.0); # plot "exact" solution
plt.ylim(yspan)
plt.show()
(0, 5)
import numpy as np
def lincvg(x):
    y = []
    for k in xrange(30):
        x = x +(cos(x)+1)/sin(x)
        show(x)
        y += [x]
    array1=np.array(y)
    err = array1 - x
    rate = err[1:]/err[:-1]
    return err , rate

(err,rate)=lincvg(3.0)
show(err,rate)
3.07091484430265\displaystyle 3.07091484430265
3.10626846715633\displaystyle 3.10626846715633
3.12393239716286\displaystyle 3.12393239716286
3.13276275488194\displaystyle 3.13276275488194
3.13717773292117\displaystyle 3.13717773292117
3.13938519684106\displaystyle 3.13938519684106
3.14048892566363\displaystyle 3.14048892566363
3.14104078968276\displaystyle 3.14104078968276
3.14131672164333\displaystyle 3.14131672164333
3.14145468761757\displaystyle 3.14145468761757
3.14152367060372\displaystyle 3.14152367060372
3.14155816209674\displaystyle 3.14155816209674
3.14157540784243\displaystyle 3.14157540784243
3.14158403071807\displaystyle 3.14158403071807
3.14158834215044\displaystyle 3.14158834215044
3.14159049786487\displaystyle 3.14159049786487
3.14159157573234\displaystyle 3.14159157573234
3.14159211464282\displaystyle 3.14159211464282
3.14159238408895\displaystyle 3.14159238408895
3.14159251879833\displaystyle 3.14159251879833
3.14159258633843\displaystyle 3.14159258633843
3.14159261935554\displaystyle 3.14159261935554
3.14159263557063\displaystyle 3.14159263557063
3.14159264173198\displaystyle 3.14159264173198
3.14159265109478\displaystyle 3.14159265109478
3.14159265109478\displaystyle 3.14159265109478
3.14159265109478\displaystyle 3.14159265109478
3.14159265109478\displaystyle 3.14159265109478
3.14159265109478\displaystyle 3.14159265109478
3.14159265109478\displaystyle 3.14159265109478
<string>:9: RuntimeWarning: invalid value encountered in divide
[ -7.06778068e-02 -3.53241839e-02 -1.76602539e-02 -8.82989621e-03 -4.41491817e-03 -2.20745425e-03 -1.10372543e-03 -5.51861412e-04 -2.75929451e-04 -1.37963477e-04 -6.89804911e-05 -3.44889980e-05 -1.72432523e-05 -8.62037671e-06 -4.30894434e-06 -2.15322991e-06 -1.07536244e-06 -5.36451957e-07 -2.67005834e-07 -1.32296451e-07 -6.47563523e-08 -3.17392375e-08 -1.55241460e-08 -9.36279854e-09 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 0.49979174 0.49994797 0.49998693 0.49999661 0.49999891 0.49999923 0.49999882 0.49999773 0.49999548 0.49999096 0.49998192 0.49996385 0.49992754 0.49985569 0.4997117 0.49941831 0.49885688 0.49772553 0.4954815 0.48947913 0.49013319 0.48911528 0.60311199 -0. nan nan nan nan nan]

cos(1.0)+1
1.54030230586814