Chapt 2: IVP of ODE
Feng Wan 
B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaoton University
Code
from sympy import *
def ode(x, y):
return x**3 + 5 * x**2.0
# Euler method for solving ODE
def euler_method(ode, x0, y0, h, num_steps):
x_values = [x0]
y_values = [y0]
for _ in range(num_steps):
x = x_values[-1]
y = y_values[-1]
y_new = y + h * ode(x, y)
x_values.append(x + h)
y_values.append(y_new)
return x_values, y_values
# Initial conditions
x0 = -3
y0 = 1
h = 1
num_steps = 5
var('x')
var('y', cls = Function)
explicit = dsolve(Eq(Derivative(y(x), x), x**3 + 5 * x**2), y(x), ics={y(x0): y0})
call_fct = lambdify((x), explicit.rhs)
# Solve ODE using Euler method
x_values, y_values = euler_method(ode, x0, y0, h, num_steps)
# Plot the solution
plt.figure()
plt.plot(x_values, y_values, marker='o', color='b', label='Euler Method')
xx = np.linspace(x0, -x0, 100)
plt.plot(xx, call_fct(xx), color = 'r', label = 'Explicit')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.legend()
plt.grid(True)
plt.show()
Exercise 1 \(dy/dx = -xy; ~~~~ y(0) = 1\), 求\(y(1), ~y(3)\)，对比不同步长\(h\); \(y = e^{-x^2/2}\)
Solution 1.
Code
from sympy import *
from sympy.printing.mathml import print_mathml
from IPython.display import Markdown
result = dsolve(Eq(Derivative(y(x), x), -x*y(x)), y(x), ics={y(0): 1})
call_fct = lambdify((x), result.rhs)
def ode(x, y):
return - x * y
def euler(x0, y0, h, xp):
steps = int((xp - x0) / h)
yi = y0
xi = x0
for _ in np.arange(steps):
xi = xi + h
yi = yi + h * ode(xi, yi)
return xi, yi
hs = [0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001]
x0 = 0
y0 = 1
xp1 = 1
xp2 = 3
print('\t%s\t%s\t%s\t%s\t%s' %("h", "y(1)_n", "y(1)", "y(3)_n", "y(3)"))
for h in hs:
xp11, yp1 = euler(x0, y0, h, xp1)
xp22, yp2 = euler(x0, y0, h, xp2)
y1 = call_fct(xp1)
y2 = call_fct(xp2)
print("%6.4f\t%6.4f\t%6.4f\t%6.4f\t%6.4f" %(h, yp1, y1, yp2, y2))
plt.figure(figsize = (12, 2))
plt.plot(xx, call_fct(xx))
plt.show()
Markdown(f'$${latex( result )}$$') h y(1)_n y(1) y(3)_n y(3)
0.5000 0.3750 0.6065 0.0000 0.0111
0.2000 0.5223 0.6065 0.0018 0.0111
0.1000 0.5653 0.6065 0.0055 0.0111
0.0500 0.5861 0.6065 0.0080 0.0111
0.0200 0.5984 0.6065 0.0098 0.0111
0.0100 0.6025 0.6065 0.0105 0.0111
0.0050 0.6045 0.6065 0.0108 0.0111
0.0020 0.6057 0.6065 0.0110 0.0111
0.0010 0.6061 0.6065 0.0110 0.0111
\[y{\left(x \right)} = e^{- \frac{x^{2}}{2}}\]
Exercise 2 1D harmonic oscillator equation: \(m\ddot{x} = -kx\) with IC: \(x_0 = 5; p_0 = 0\)
✅ Explicit solution: \[ x = A\cos(\omega_0 t + \phi); ~~~ p = -mA\omega_0\sin(\omega_0t+\phi) \] with \(\omega_0 \equiv \sqrt{\frac{k}{m}}\).
\(A\) & \(\phi\) are determined by IC: \(\frac{1}{2}kx_0^2 + \frac{1}{2}mv_0^2 = \frac{1}{2}kA^2;~~~\phi=\cos^{-1}(x_0/A)\).
\[p{\left(t \right)} = 0.790569415042095 i e^{0.316227766016838 i t} - 0.790569415042095 i e^{- 0.316227766016838 i t}\]
Code
\[\left[ y{\left(t \right)} = \frac{\sin{\left(2 \sqrt{\pi} t \right)}}{2 \sqrt{\pi}}, \ p{\left(t \right)} = \cos{\left(2 \sqrt{\pi} t \right)}\right]\]
Code
def ode(x, y):
res = y.copy()
y1 = res[0]
y2 = res[1]
z1 = res[2]
z2 = res[3]
y11 = z1
y22 = z2
z11 = -(z1 + 0.2 * (z1 - z2)) + np.sin(y1) + np.sin(x)
z22 = -(z2 + 0.2 * (z2 - z1)) + np.cos(y2) + np.sin(x)
return np.array([y11, y22, z11, z22])
def rk4(func, x0, y0, h, xp):
steps = int((xp - x0) / h)
result = np.zeros((steps + 1, len(y0)))
xis = np.zeros((steps + 1, 1))
xi = x0
yi = y0
xis[0] = x0
result[0, :] = yi
for i in range(steps):
k1 = func(xi, yi)
k2 = func(xi + 0.5 * h, yi + 0.5 * h * k1)
k3 = func(xi + 0.5 * h, yi + 0.5 * h * k2)
k4 = func(xi + h, yi + h * k3)
xi = xi + h
yi = yi + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
xis[i + 1] = xi
result[i + 1, :] = yi
return xis, result
h = 0.01
x0 = 0
xp = 30.0
y0 = np.array([1.0, 0.0, 2.0, 5.0])
xi, yi = rk4(ode, x0, y0, h, xp)
plt.plot(xi, yi[:, 0], label = r'$y_1$')
plt.plot(xi, yi[:, 1], label = r'$y_2$')
plt.plot(xi, yi[:, 2], label = r'$dy_1/dt$')
plt.plot(xi, yi[:, 3], label = r'$dy_2/dt$')
plt.xlabel(r'$t$')
plt.legend()
plt.show()
Code
def func(y, t):
return [y[1], -4*np.pi**2*y[0]]
y0 = [0, 1]
t = np.linspace(0, 10, 1001)
sol = scipy.integrate.odeint(func, y0, t)
plt.figure()
plt.subplot(121)
plt.plot(t, sol[:, 0], 'r', label='y')
plt.plot(t, sol[:, 1], 'b', label='p')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.subplot(122)
plt.plot(sol[:, 0], sol[:, 1], 'b', label='y-p')
plt.legend(loc='best')
plt.xlabel('y')
plt.ylabel('p')
plt.grid()
plt.show()
Code
def func(y, t):
return [y[1], y[2], y[2] + t]
y0 = [8, 7, 4]
x0 = 1
t = np.linspace(1, 10, 1001)
sol = scipy.integrate.odeint(func, y0, t)
plt.figure()
plt.plot(t, sol[:, 0], 'r', label='y')
callfct = lambdify((x), res.rhs)
plt.plot(t, callfct(t), 'b-.', label = 'explicit')
plt.legend(loc='best')
plt.xlabel('x')
plt.grid()
plt.show()
Code
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
# Parameters
x_res, y_res = 300, 300
xmin, xmax = -1.5, 1.5
width = xmax - xmin
ymin, ymax = -1.5, 1.5
height = ymax - ymin
z_abs_max = 10
max_iter = 1000
def julia_set(c):
# Initialise an empty array (corresponding to pixels)
julia = np.zeros((x_res, y_res))
# Loop over each pixel
for ix in range(x_res):
for iy in range(y_res):
# Map pixel position to a point in the complex plane
z = complex(ix / x_res * width + xmin,
iy / y_res * height + ymin)
# Iterate
iteration = 0
while abs(z) <= z_abs_max and iteration < max_iter:
z = z**2 + c
iteration += 1
iteration_ratio = iteration / max_iter
# Set the pixel value to be equal to the iteration_ratio
julia[ix, iy] = iteration_ratio
# Plot the array using matplotlib's imshow
fig, ax = plt.subplots()
ax.imshow(julia, interpolation='nearest', cmap=cm.gnuplot2)
plt.axis('off')
plt.show()
c = -0.1-0.65j
julia_set(c)
一根长度为 \(l\) 的钟摆被限制在一垂直的平面内，在强迫外力 \(f_d\) 和阻力 \(f_r\) 的作用下振荡运动。钟摆的运动可以通过 Newton 方程来描述 \[ ma = f_g + f_d + f_r \] 其中 \(f_g = -mg \sin \theta\) 是重力在运动方向的分力，\(a=l d^2\theta/dt^2\) 是沿切线方向的加速度， \(\theta\) 是杆和垂线的夹角。
Code
q = 0.5
b = 1.9
omega = 2/3.0
theta0 = np.pi * 160 / 180
y0 = np.array([theta0, 0])
t_end = 1000
steps = 160000
h = t_end / steps
def ode(x, y):
x1 = y[0]
x2 = y[1]
res = np.array([x2, -q * x2 - np.sin(x1) + b * np.cos(omega * x)])
return res
xi, yi = rk4(ode, 0, y0, h, t_end)
plt.plot(xi, yi[:, 0])
plt.plot(xi, yi[:, 1])
plt.show()
Code
q = 0
b = 0
omega = 2/3.0
theta0 = np.pi * 3 / 180
y0 = np.array([theta0, 0])
t_end = 50
steps = 20000
h = t_end / steps
xi, yi = rk4(ode, 0, y0, h, t_end)
plt.subplot(121)
plt.plot(xi, yi[:, 0], label = r'$\theta$')
plt.plot(xi, yi[:, 1], label = r'$\dot{\theta}$')
plt.xlabel(r'$t$')
plt.legend()
plt.subplot(122)
plt.plot(yi[:, 0], yi[:, 1], label = r'$\theta-\dot{\theta}$')
plt.xlabel(r'$\theta$')
#plt.ylabel(r'$\dot{\theta}$')
plt.legend()
plt.show()
Code
\(q=0,b=0, t=50\)，初始角度小于\(5^\circ\)。即无阻尼，小角度释放的情况，钟摆运动是一个简谐运动, i.e., \(\theta(t) = \theta_0 \cos \omega t\).
Code
q = 0
b = 0
omega = 2/3.0
theta0 = np.pi * 160 / 180
y0 = np.array([theta0, 0])
t_end = 50
steps = 20000
h = t_end / steps
xi, yi = rk4(ode, 0, y0, h, t_end)
plt.subplot(121)
plt.plot(xi, yi[:, 0], label = r'$\theta$')
plt.plot(xi, yi[:, 1], label = r'$\dot{\theta}$')
plt.xlabel(r'$t$')
plt.legend()
plt.subplot(122)
plt.plot(yi[:, 0], yi[:, 1], label = r'$\theta-\dot{\theta}$')
plt.xlabel(r'$\theta$')
#plt.ylabel(r'$\dot{\theta}$')
plt.legend()
plt.show()
Code
q = 0
b = 0
omega = 2/3.0
t_end = 50
steps = 20000
h = t_end / steps
theta0s = np.linspace(0.01, 0.99, 100) * np.pi
dur = np.zeros(theta0s.shape)
for i in range(len(theta0s)):
theta0 = theta0s[i]
y0 = np.array([theta0, 0])
xi, yi = rk4(ode, 0, y0, h, t_end)
mytheta = yi[:, 0]
mythetam = mytheta[0:-1]
mythetap = mytheta[1:]
indexes = np.where(mythetam < mythetap)
#print(i, indexes)
dur[i] = 2 * xi[indexes[0][0]]
plt.plot(theta0s, dur, '-o')
plt.xlabel(r'$\theta_0$')
plt.ylabel(r'$T$')
plt.show()
Code
q = 0
b = 0
omega = 2/3.0
theta0 = np.pi * 180 / 180
y0 = np.array([theta0, 0])
t_end = 20
steps = 80000
h = t_end / steps
xi, yi = rk4(ode, 0, y0, h, t_end)
plt.subplot(311)
plt.plot(xi, yi[:, 0], label = r'$\theta$')
plt.xlabel(r'$t$')
plt.ylabel(r'$\theta$')
plt.subplot(312)
plt.plot(xi, yi[:, 1], label = r'$\dot{\theta}$')
plt.xlabel(r'$t$')
plt.ylabel(r'$\dot{\theta}$')
plt.subplot(313)
plt.plot(yi[:, 0], yi[:, 1], label = r'$\theta-\dot{\theta}$')
plt.xlabel(r'$\theta$')
plt.ylabel(r'$\dot{\theta}$')
plt.show()
单摆做旋转运动，由初速度决定运动状态
