ODE of Joy

30 days ago by solin

ODE of Joy

This worksheet is a bridge between a basic Python tutorial and elementary numerical methods for Ordinary Differential Equations (ODE). We will learn how to loop over arrays of points, create pointwise defined functions, and plot graphs. Then we will implement the explicit Euler method and study its properties.

Creating arrays of integers:
list = range(5) print list 
       
list = range(10) for i in list: print i*i 
       

Creating arrays of real numbers:
from pylab import arange x_min = 0. x_max = 5. n = 10 h = (x_max - x_min)/n x_array = arange(x_min, x_max + h/2., h) print x_array 
       

Creating pointwise defined functions:
def f(x): return 1/(1 + x*x) y_array = [f(x) for x in x_array] print y_array 
       

Plotting pointwise defined functions:
from pylab import plot, savefig, clf, grid, legend, show clf() label = "This is my function" plot(x_array, y_array, "r-o", label=label) legend() #grid(True) savefig("a.png") 
       

Solving an ODE of the form Y' = F(Y,t), Y(0) = y0 in the interval (0, T) with the explicit Euler method:
# Importing libraries: from pylab import arange, zeros, plot, savefig, clf, grid, legend, array from pylab import pi, exp, sin, cos # Explicit Euler (RK1) method def RK1(F, y0, t_array, n): y_array = zeros(n+1) y_array[0] = y0 for i in range(n): y_array[i+1] = y_array[i] + dt*F(y_array[i], t_array[i]) return y_array # RK2 method def RK2(F, y0, t_array, n): b1 = 1./4 b2 = 3./4 c2 = 2./3 a21 = 2./3 y_array = zeros(n+1) y_array[0] = y0 for i in range(n): k1 = F(y_array[i], t_array[i]) k2 = F(y_array[i] + a21*dt*k1, t_array[i] + c2*dt) y_array[i+1] = y_array[i] + dt*(b1*k1 + b2*k2) return y_array # RK4 method def RK4(F, y0, t_array, n): b1 = 1./6 b2 = 1./3 b3 = 1./3 b4 = 1./6 c2 = 1./2 c3 = 1./2 c4 = 1. a21 = 1./2 a31 = 0. a32 = 1./2 a41 = 0. a42 = 0. a43 = 1. y_array = zeros(n+1) y_array[0] = y0 for i in range(n): k1 = F(y_array[i], t_array[i]) k2 = F(y_array[i] + dt*(a21*k1), t_array[i] + c2*dt) k3 = F(y_array[i] + dt*(a31*k1 + a32*k2), t_array[i] + c3*dt) k4 = F(y_array[i] + dt*(a41*k1 + a42*k2 + a43*k3), t_array[i] + c4*dt) y_array[i+1] = y_array[i] + dt*(b1*k1 + b2*k2 + b3*k3 + b4*k4) return y_array # Plotting exact solution def plot_exact(E, t_min, t_max, n_exact): dt_exact = (t_max - t_min)/n_exact t_array_exact = arange(t_min, t_max + dt_exact/2., dt_exact) E_array = array([E(t) for t in t_array_exact]) label = "Exact solution E(t)" plot(t_array_exact, E_array, "g-", label=label) # Main program: # Problem definition: t_min = 0. t_max = pi K = 3 # Right-hand side def F(y,t): return K*cos(K*t) #return -y # Initial condition y0 = 0. # Is exact solution available? exact_known = True # Exact solution: def E(t): #return 3*exp(-t) return sin(K*t) # Division n = 30 # Time step dt = (t_max - t_min)/n # Array of discrete time points t_array = arange(t_min, t_max + dt/2., dt) clf() # Explicit Euler (RK1) method Y_RK1_array = RK1(F, y0, t_array, n) label = "RK1" plot(t_array, Y_RK1_array, "b-", label=label) # Second-order Runge-Kutta (RK2) method Y_RK2_array = RK2(F, y0, t_array, n) label = "RK2" plot(t_array, Y_RK2_array, "y-", label=label) # Fourth-order Runge-Kutta (RK4) method Y_RK4_array = RK4(F, y0, t_array, n) label = "RK4" plot(t_array, Y_RK4_array, "r-", label=label) # Plot exact solution if known if exact_known == True: plot_exact(E, t_min, t_max, 500) legend() grid(True) savefig("a.png")