• API
• FAQ
• Tools
• Archive
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.

Top