SHARE
TWEET

Untitled

a guest Dec 17th, 2018 68 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/bin/env python
  2. """
  3. Find the solution for the second order differential equation
  4. u'' = -u
  5. with u(0) = 10 and u'(0) = -5
  6. using the Euler and the Runge-Kutta methods.
  7. This works by splitting the problem into 2 first order differential equations
  8. u' = v
  9. v' = f(t,u)
  10. with u(0) = 10 and v(0) = -5
  11. """
  12.  
  13. from math import cos, sin
  14.  
  15. def f(t, u):
  16.     return -u
  17.  
  18. def exact(u0, du0, t):
  19.     # analytical solution
  20.     return u0 * cos(t) + du0 * sin(t)
  21.  
  22. def iterate(func, u, v, tmax, n):
  23.     dt = tmax/(n-1)
  24.     t = 0.0
  25.  
  26.     for i in range(n):
  27.         u,v = func(u,v,t,dt)
  28.         t += dt
  29.  
  30.     return u
  31.  
  32. def euler_iter(u, v, t, dt):
  33.     v_new = v + dt * f(t, u)
  34.     u_new = u + dt * v
  35.     return u_new, v_new
  36.  
  37. def rk_iter(u, v, t, dt):
  38.     k1 = f(t,u)
  39.     k2 = f(t+dt*0.5,u+k1*0.5*dt)
  40.     k3 = f(t+dt*0.5,u+k2*0.5*dt)
  41.     k4 = f(t+dt,u+k3*dt)
  42.  
  43.     v += dt * (k1+2*k2+2*k3+k4)/6
  44.  
  45.     # v doesn't explicitly depend on other variables
  46.     k1 = k2 = k3 = k4 = v
  47.  
  48.     u += dt * (k1+2*k2+2*k3+k4)/6
  49.  
  50.     return u,v
  51.  
  52. euler = lambda u, v, tmax, n: iterate(euler_iter, u, v, tmax, n)
  53. runge_kutta = lambda u, v, tmax, n: iterate(rk_iter, u, v, tmax, n)
  54.  
  55. def plot_result(u, v, tmax, n):
  56.     dt = tmax/(n-1)
  57.     t = 0.0
  58.     allt = []
  59.     error_euler = []
  60.     error_rk = []
  61.     r_exact = []
  62.     r_euler = []
  63.     r_rk = []
  64.  
  65.     u0 = u_euler = u_rk = u
  66.     v0 = v_euler = v_rk = v
  67.  
  68.     for i in range(n):
  69.         u = exact(u0, v0, t)
  70.         u_euler, v_euler = euler_iter(u_euler, v_euler, t, dt)
  71.         u_rk, v_rk = rk_iter(u_rk, v_rk, t, dt)
  72.         allt.append(t)
  73.         error_euler.append(abs(u_euler-u))
  74.         error_rk.append(abs(u_rk-u))
  75.         r_exact.append(u)
  76.         r_euler.append(u_euler)
  77.         r_rk.append(u_rk)
  78.         t += dt
  79.  
  80.     _plot("error.png", "Error", "time t", "error e", allt, error_euler, error_rk)
  81.     #_plot("result.png", "Result", "time t", "u(t)", allt, r_euler, r_rk, r_exact)
  82.  
  83. def _plot(out, title, xlabel, ylabel, allt, euler, rk, exact=None):
  84.     import matplotlib.pyplot as plt
  85.  
  86.     plt.title(title)
  87.  
  88.     plt.ylabel(ylabel)
  89.     plt.xlabel(xlabel)
  90.  
  91.     plt.plot(allt, euler, 'b-', label="Euler")
  92.     plt.plot(allt, rk, 'r--', label="Runge-Kutta")
  93.  
  94.     if exact:
  95.         plt.plot(allt, exact, 'g.', label='Exact')
  96.  
  97.     plt.legend(loc=4)
  98.     plt.grid(True)
  99.  
  100.     plt.savefig(out, dpi=None, facecolor='w', edgecolor='w',
  101.                 orientation='portrait', papertype=None, format=None,
  102.                 transparent=False)
  103.  
  104. u0 = 10
  105. du0 = v0 = -5
  106. tmax = 10.0
  107. n = 2000
  108.  
  109. print "t=", tmax
  110. print "euler =", euler(u0, v0, tmax, n)
  111. print "runge_kutta=", runge_kutta(u0, v0, tmax, n)
  112. print "exact=", exact(u0, v0, tmax)
  113.  
  114. plot_result(u0, v0, tmax*2, n*2)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top