25/10/2017

Lab 11: Numerical Solutions of ODEs

---

In this lab, you will implement the different Runge-Kutta Methods we studied in class. Please test your implementation on the function given in this webpage. Compare your answers with the ones given below.

from math import *
from matplotlib import pyplot as plt
import numpy as np
import ode_solvers as ode

np.set_printoptions(precision=5, suppress=True)

def ode_euler (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

def ode_rk1 (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

def ode_rk2_heun (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

def ode_rk2_midpoint (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

def ode_rk2_ralston (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

# This is THE Runge-Kutta Method, rk4
def ode_solve (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

def ode_rk4_38 (f, tn, tdelta, x0=0, t0=0) :
    return tn, xn

# Use your code to solve the following differential eqn
#         dx/dt + 2x = cos(4t)  x(0) = 3
# Compare the solutions for different methods. 
# The exact solution is
#         x(t) = 0.2sin(4t)+0.1cos(4t)+2.9exp(-2t)

xt = np.zeros((51,2))
count = 0
for i in np.linspace(0.0, 10.0, 51) :
    xt[count,:] = ode.ode_rk4_16(ode.xdot, i, 0.2, x0=3.)
    count = count + 1

print xt[:,1]

[ 3.       2.15726  1.50029  0.93507  0.47426  0.17601  0.07285  0.12785
  0.24096  0.29886  0.23644  0.07148 -0.10944 -0.2056  -0.16473 -0.01568
  0.14842  0.2262   0.16926  0.01132 -0.15237 -0.22288 -0.15769  0.00349
  0.16278  0.22348  0.14873 -0.01618 -0.17122 -0.22238 -0.13862  0.02924
  0.17937  0.2207   0.12817 -0.04211 -0.18684 -0.21823 -0.11725  0.05486
  0.19369  0.21503  0.10594 -0.06741 -0.19988 -0.2111  -0.09427  0.07974
  0.20538  0.20644  0.08227]

def exact (x) :
    return 0.2 * sin(4.0 * x) + 0.1 * cos(4.0 * x) + 2.9 * exp(-2.0 * x)

count = 0
xte = np.zeros(51)
for i in np.linspace(0.0, 10.0, 51) :
    xte[count] = exact(i)
    count = count + 1
    
print xte

[ 3.       2.15707  1.50005  0.93482  0.474    0.17575  0.0726   0.12765
  0.24084  0.29881  0.23644  0.07148 -0.10947 -0.20567 -0.16481 -0.01574
  0.1484   0.22623  0.16931  0.01136 -0.15237 -0.22293 -0.15776  0.00344
  0.16278  0.22353  0.14879 -0.01614 -0.17123 -0.22242 -0.13868  0.0292
  0.17938  0.22075  0.12823 -0.04208 -0.18686 -0.21829 -0.11731  0.05483
  0.19371  0.21509  0.106   -0.06739 -0.1999  -0.21115 -0.09432  0.07972
  0.20541  0.2065   0.08233]

plt.plot(xt[:,0], xt[:,1], 'g-')
plt.plot(xt[:,0], xte, 'b.')
plt.show()