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 odeint
from scipy.interpolate import interp1d
from sage.plot.colors import rainbow
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):
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)

$\displaystyle 3.07091484430265$
$\displaystyle 3.10626846715633$
$\displaystyle 3.12393239716286$
$\displaystyle 3.13276275488194$
$\displaystyle 3.13717773292117$
$\displaystyle 3.13938519684106$
$\displaystyle 3.14048892566363$
$\displaystyle 3.14104078968276$
$\displaystyle 3.14131672164333$
$\displaystyle 3.14145468761757$
$\displaystyle 3.14152367060372$
$\displaystyle 3.14155816209674$
$\displaystyle 3.14157540784243$
$\displaystyle 3.14158403071807$
$\displaystyle 3.14158834215044$
$\displaystyle 3.14159049786487$
$\displaystyle 3.14159157573234$
$\displaystyle 3.14159211464282$
$\displaystyle 3.14159238408895$
$\displaystyle 3.14159251879833$
$\displaystyle 3.14159258633843$
$\displaystyle 3.14159261935554$
$\displaystyle 3.14159263557063$
$\displaystyle 3.14159264173198$
$\displaystyle 3.14159265109478$
$\displaystyle 3.14159265109478$
$\displaystyle 3.14159265109478$
$\displaystyle 3.14159265109478$
$\displaystyle 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