| Hosted by CoCalc | Download
# -*- 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 range(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
/home/user/cloud-examples/sage/2018-02-15-162612_Picard_Numérico_Integral_y_ODE.sagews:9: RuntimeWarning: invalid value encountered in true_divide from sage.plot.colors import rainbow
[-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