taylorClass.py
#
# Class and explicit integrators for ODE
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
from matplotlib import pyplot as plt
# ============================= a so simple class ! =======================
class ExplMethods(object):
"Class of explicit integrators for ODE equations"
#
# Un constructeur qui est exécuté automatiquement lorsqu'on instancie un
# nouvel objet de la classe : il s'appelle obligatoirement __init__()
#
def __init__(self,name,color,f):
self.name = name
self.f = f
self.color = color
# ============================= mainProgram ===============================
integrators = [ExplMethods("Explicit Euler order 1",".-b",
lambda u,x,h : h*(x-u)/2 ),
ExplMethods("Explicit Taylor order 4",".-r",
lambda u,x,h : h*(x-u)/2+(2-x+u)*((h**2)/8-(h**3)/48+(h**4)/384) ) ]
u = lambda x : 3*exp(-x/2)-2+x
Xstart = 0; Xend = 3; Ustart = 1;
for integrator in integrators:
print(" ====== %s method=====" % integrator.name)
error = zeros(4)
for j in range(4):
n = 3*pow(2,j)
h = (Xend-Xstart)/n
X = linspace(Xstart,Xend,n+1)
U = zeros(n+1); U[0] = Ustart
for i in range(n):
U[i+1] = U[i] + integrator.f(U[i],X[i],h)
plt.plot(X,U,integrator.color)
error[j] = abs(U[-1]-u(Xend))
print(" ==== %s : h=%5.3f : eh(Xend) = %8.2e " % (integrator.name,h,error[j]))
order = mean(log(error[:-1]/error[1:])/log(2))
print(" ============= Estimated order : %.4f " % order)
plt.show()