应用修改后的欧拉求解 Python 中的钟摆 ODE
Applying Modified Euler to solve a Pendulum ODE in Python
所以我试图在 python 中实现一些数值方法，但我遇到了一些问题，我的所有函数输出的内容与常规欧拉方法大致相同。我认为这是因为我在将方法实现到代码中时以某种方式搞砸了。
我的钟摆是这样定义的：
def func(y,t):
### Simplified the Function to remove friction since it cancelled out
x,v = y[:3],y[3:6]
grav = np.array([0., 0., -9.8 ])
lambd = (grav.dot(x)+v.dot(v))/x.dot(x)
return np.concatenate([v, grav - lambd*x] )
def dF_matrix(y):
n=y.size
dF=np.zeros((6,6))
xp=np.array([y[1],y[2],y[3]])[np.newaxis]
mass=1.
F1=0.
F2=0.
F3=-mass*9.8
F=np.array([F1,F2,F3])[np.newaxis]
phix=2.*y[0]
phiy=2.*y[4]
phiz=2.*y[5]
G=np.array([phix,phiy,phiz])[np.newaxis]
H=2.*np.eye(3)
lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)
dF[0,3]=1
dF[1,4]=1
dF[2,5]=1
dF[3,0]=(y[0]*F1+2*lambd)/mass
dF[3,1]=(y[0]*F2)/mass
dF[3,2]=(y[0]*F3)/mass
dF[3,3]=phix*y[1]
dF[3,4]=phix*y[2]
dF[3,5]=phix*y[3]
dF[4,0]=(y[4]*F1)/mass
dF[4,1]=(y[4]*F2+2*lambd)/mass
dF[4,2]=(y[4]*F3)/mass
dF[4,3]=phiy*y[1]
dF[4,4]=phiy*y[2]
dF[4,5]=phiy*y[3]
dF[5,0]=(y[5]*F1)/mass
dF[5,1]=(y[5]*F2)/mass
dF[5,2]=(y[5]*F3+2*lambd)/mass
dF[5,3]=phiz*y[1]
dF[5,4]=phiz*y[2]
dF[5,5]=phiz*y[3]
return dF
为了整合ODE函数，我做了如下的功能（在上一个帖子中得到了其他人的帮助）：
from scipy.integrate import odeint
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
正向欧拉法
def forward_Euler(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
y[i + 1, :] = y[i, :] + np.asarray(function(y[i, :], time[i])) * dt
return y
修正欧拉法
错误从这里开始
我得到的错误是：
RuntimeWarning：double_scalars 中遇到无效值
lambd = (grav.dot(x)+v.dot(v))/x.dot(x)
def modified_Euler(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix))) # creates the matrix that we will fill
y[0, :] = y_matrix # sets the initial values of the matrix
for i in range(len(time) - 1): # apply the Euler
dt = time[i + 1] - time[i] # Step size
k1 = np.asarray(function(y[i, :], time[i])*dt)
k2 = np.asarray(function(y[i] + k1, time[i+1])*dt)
y[i + 1, :] = y[i, :] + .5 * (k1 + k2)
return y
Adams-Bashforth 二阶
def Adams_Bash_2nd(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
dt = time[1] - time[0]
f_0 = function(y[0], time[0])
y[1] = y[0] + dt * f_0
y[1] = y[0] + 0.5*dt * (f_0+function(y[1], time[1]))
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
f_1 = function(y[i, :], time[i])
f_2 = function(f_1-1, time[i-1])
y[i + 1] = y[i] + 0.5 * dt * (3 * f_1 - f_2)
return y
Adams Bashforth Moulton 方法
def Adams_Moulton(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
### predictor formula
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
f_1 = function(y[i, :], time[i])
f_2 = function(f_1-1, time[i-1])
y[i + 1, :] = y[i, :] + dt * f_1 + ((dt**2)/2) * f_2
### Corrector formula
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
k_1 = 9 * (function(y[i, :], time[i+1]))
k_2 = 19 * (function(y[i, :], time[i]))
k_3 = 5 * (function(y[i, :], time[i-1]))
k_4 = (function(y[i, :], time[i-2]))
y[i + 1, :] = y[i] + (dt/24) * (k_1 + k_2 - k_3 + k_4)
return y
在下一个函数中使用的 RK4 步骤
def RK4_step(f,y,t,dt, N=1):
dt /= N;
for k in range(N):
k1=f(y,t)*dt; k2=f(y+k1/2,t+dt/2)*dt; k3=f(y+k2/2,t+dt/2)*dt; k4=f(y+k3,t+dt)*dt;
y, t = y+(k1+2*(k2+k3)+k4)/6, t+dt
return y
Adams-Bashforth Moulton 方法四阶
def Adams_Moulton_4th(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0] = y_matrix
### bootstrap steps with 4th order one-step method
dt = time[4] - time[0]
y[4] = RK4_step(function, y[0], time[0], dt, N=4)
y[5] = RK4_step(function, y[4], time[4], dt, N=4)
y[1] = RK4_step(function, y[5], time[5], dt, N=4)
f_m2 = function(y[0], time[0])
f_m1 = function(y[4], time[4])
f_0 = function(y[5], time[5])
f_1 = function(y[1], time[1])
for i in range(3, len(time) - 1):
### predictor formula 4th order [ 55/24, -59/24, 37/24, -3/8 ]
f_m3, f_m2, f_m1, f_0 = f_m2, f_m1, f_0, f_1
y[i + 1] = y[i] + (dt / 24) * (55 * f_0 - 59 * f_m1 + 37 * f_m2 - 9 * f_m3)
f_1 = function(y[i + 1], time[i + 1])
### Corrector formula 4th order [ 3/8, 19/24, -5/24, 1/24 ]
y[i + 1] = y[i] + (dt / 24) * (9 * f_1 + 19 * f_0 - 5 * f_m1 + f_m2)
f_1 = function(y[i + 1], time[i + 1])
return y
我决定将我测试函数的方式编程为一个函数，该函数从上一次迭代中删除了大量行
# initial condition
y0 = np.array([0.0, 1.0, 0.0, 0.8, 0.0, 1.2])
def test_function(test_function):
print(test_function.__name__ + "...")
nt = 2500
time_start = process_time()
# time points
t = np.linspace(0, 25, nt)
# solve ODE
y1 = test_function(func, y0, t)
time_elapsed = (process_time() - time_start)
print('elapsed time', time_elapsed)
# compute residual:
r = y1[:, 0] ** 2 + y1[:, 1] ** 2 + y1[:, 2] ** 2 - 1
rmax1 = np.max(np.abs(r))
print('error', rmax1)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(y1[:, 0], y1[:, 1], y1[:, 2], 'gray')
plt.show()
test_function(odeint)
test_function(forward_Euler)
test_function(modified_Euler)
test_function(Adams_Bash_2nd)
test_function(Adams_Moulton)
test_function(Adams_Moulton_4th)
mod化欧拉方法没有步骤i -> i+1之外的访问点，没有i-1（注意你的源文档 python 代码而非公式中的步骤是 i-1 -> i，循环从适当增加的索引开始）。它只是（因为你到处都可以找到 mod。讨论了 Euler 或 Heun 方法）
k1 = f(y[i] , t[i ])*dt;
k2 = f(y[i]+k1, t[i+1])*dt;
y[i+1] = y[i] + 0.5*(k1+k2);
相比之下，2 阶的 Adams-Bashford 方法和大于 2 阶的 Adams-Moulton 方法 Do 从步骤 i -> i+1 之前访问点，形式上是一个在 AB2
y[i+1] = y[i] + 0.5*dt * (3*f[i] - f[i-1])
对于第一个实现，以与 y 数组相同的方式声明 f 数组以逐字实现此公式是有意义的。仅保留一个 f 值的短数组可能更经济，这些值被移动或旋转以访问最后几个 f 值。
注意您需要使用其他类似或更高阶的方法初始化y[1]和f[1]。或者如果你想有一个 "pure" 运行 的方法，你需要初始化 y[-1] 和 f[-1] 并进一步返回以便 y[1] 可以用方法公式。
所以我试图在 python 中实现一些数值方法，但我遇到了一些问题，我的所有函数输出的内容与常规欧拉方法大致相同。我认为这是因为我在将方法实现到代码中时以某种方式搞砸了。
我的钟摆是这样定义的：
def func(y,t):
### Simplified the Function to remove friction since it cancelled out
x,v = y[:3],y[3:6]
grav = np.array([0., 0., -9.8 ])
lambd = (grav.dot(x)+v.dot(v))/x.dot(x)
return np.concatenate([v, grav - lambd*x] )
def dF_matrix(y):
n=y.size
dF=np.zeros((6,6))
xp=np.array([y[1],y[2],y[3]])[np.newaxis]
mass=1.
F1=0.
F2=0.
F3=-mass*9.8
F=np.array([F1,F2,F3])[np.newaxis]
phix=2.*y[0]
phiy=2.*y[4]
phiz=2.*y[5]
G=np.array([phix,phiy,phiz])[np.newaxis]
H=2.*np.eye(3)
lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)
dF[0,3]=1
dF[1,4]=1
dF[2,5]=1
dF[3,0]=(y[0]*F1+2*lambd)/mass
dF[3,1]=(y[0]*F2)/mass
dF[3,2]=(y[0]*F3)/mass
dF[3,3]=phix*y[1]
dF[3,4]=phix*y[2]
dF[3,5]=phix*y[3]
dF[4,0]=(y[4]*F1)/mass
dF[4,1]=(y[4]*F2+2*lambd)/mass
dF[4,2]=(y[4]*F3)/mass
dF[4,3]=phiy*y[1]
dF[4,4]=phiy*y[2]
dF[4,5]=phiy*y[3]
dF[5,0]=(y[5]*F1)/mass
dF[5,1]=(y[5]*F2)/mass
dF[5,2]=(y[5]*F3+2*lambd)/mass
dF[5,3]=phiz*y[1]
dF[5,4]=phiz*y[2]
dF[5,5]=phiz*y[3]
return dF
为了整合ODE函数，我做了如下的功能（在上一个帖子中得到了其他人的帮助）：
from scipy.integrate import odeint
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
正向欧拉法
def forward_Euler(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
y[i + 1, :] = y[i, :] + np.asarray(function(y[i, :], time[i])) * dt
return y
修正欧拉法
错误从这里开始
我得到的错误是： RuntimeWarning：double_scalars 中遇到无效值 lambd = (grav.dot(x)+v.dot(v))/x.dot(x)
def modified_Euler(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix))) # creates the matrix that we will fill
y[0, :] = y_matrix # sets the initial values of the matrix
for i in range(len(time) - 1): # apply the Euler
dt = time[i + 1] - time[i] # Step size
k1 = np.asarray(function(y[i, :], time[i])*dt)
k2 = np.asarray(function(y[i] + k1, time[i+1])*dt)
y[i + 1, :] = y[i, :] + .5 * (k1 + k2)
return y
Adams-Bashforth 二阶
def Adams_Bash_2nd(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
dt = time[1] - time[0]
f_0 = function(y[0], time[0])
y[1] = y[0] + dt * f_0
y[1] = y[0] + 0.5*dt * (f_0+function(y[1], time[1]))
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
f_1 = function(y[i, :], time[i])
f_2 = function(f_1-1, time[i-1])
y[i + 1] = y[i] + 0.5 * dt * (3 * f_1 - f_2)
return y
Adams Bashforth Moulton 方法
def Adams_Moulton(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0, :] = y_matrix
### predictor formula
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
f_1 = function(y[i, :], time[i])
f_2 = function(f_1-1, time[i-1])
y[i + 1, :] = y[i, :] + dt * f_1 + ((dt**2)/2) * f_2
### Corrector formula
for i in range(len(time) - 1):
dt = time[i + 1] - time[i]
k_1 = 9 * (function(y[i, :], time[i+1]))
k_2 = 19 * (function(y[i, :], time[i]))
k_3 = 5 * (function(y[i, :], time[i-1]))
k_4 = (function(y[i, :], time[i-2]))
y[i + 1, :] = y[i] + (dt/24) * (k_1 + k_2 - k_3 + k_4)
return y
在下一个函数中使用的 RK4 步骤
def RK4_step(f,y,t,dt, N=1):
dt /= N;
for k in range(N):
k1=f(y,t)*dt; k2=f(y+k1/2,t+dt/2)*dt; k3=f(y+k2/2,t+dt/2)*dt; k4=f(y+k3,t+dt)*dt;
y, t = y+(k1+2*(k2+k3)+k4)/6, t+dt
return y
Adams-Bashforth Moulton 方法四阶
def Adams_Moulton_4th(function, y_matrix, time):
y = np.zeros((np.size(time), np.size(y_matrix)))
y[0] = y_matrix
### bootstrap steps with 4th order one-step method
dt = time[4] - time[0]
y[4] = RK4_step(function, y[0], time[0], dt, N=4)
y[5] = RK4_step(function, y[4], time[4], dt, N=4)
y[1] = RK4_step(function, y[5], time[5], dt, N=4)
f_m2 = function(y[0], time[0])
f_m1 = function(y[4], time[4])
f_0 = function(y[5], time[5])
f_1 = function(y[1], time[1])
for i in range(3, len(time) - 1):
### predictor formula 4th order [ 55/24, -59/24, 37/24, -3/8 ]
f_m3, f_m2, f_m1, f_0 = f_m2, f_m1, f_0, f_1
y[i + 1] = y[i] + (dt / 24) * (55 * f_0 - 59 * f_m1 + 37 * f_m2 - 9 * f_m3)
f_1 = function(y[i + 1], time[i + 1])
### Corrector formula 4th order [ 3/8, 19/24, -5/24, 1/24 ]
y[i + 1] = y[i] + (dt / 24) * (9 * f_1 + 19 * f_0 - 5 * f_m1 + f_m2)
f_1 = function(y[i + 1], time[i + 1])
return y
我决定将我测试函数的方式编程为一个函数，该函数从上一次迭代中删除了大量行
# initial condition
y0 = np.array([0.0, 1.0, 0.0, 0.8, 0.0, 1.2])
def test_function(test_function):
print(test_function.__name__ + "...")
nt = 2500
time_start = process_time()
# time points
t = np.linspace(0, 25, nt)
# solve ODE
y1 = test_function(func, y0, t)
time_elapsed = (process_time() - time_start)
print('elapsed time', time_elapsed)
# compute residual:
r = y1[:, 0] ** 2 + y1[:, 1] ** 2 + y1[:, 2] ** 2 - 1
rmax1 = np.max(np.abs(r))
print('error', rmax1)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(y1[:, 0], y1[:, 1], y1[:, 2], 'gray')
plt.show()
test_function(odeint)
test_function(forward_Euler)
test_function(modified_Euler)
test_function(Adams_Bash_2nd)
test_function(Adams_Moulton)
test_function(Adams_Moulton_4th)
mod化欧拉方法没有步骤i -> i+1之外的访问点，没有i-1（注意你的源文档 python 代码而非公式中的步骤是 i-1 -> i，循环从适当增加的索引开始）。它只是（因为你到处都可以找到 mod。讨论了 Euler 或 Heun 方法）
k1 = f(y[i] , t[i ])*dt;
k2 = f(y[i]+k1, t[i+1])*dt;
y[i+1] = y[i] + 0.5*(k1+k2);
相比之下，2 阶的 Adams-Bashford 方法和大于 2 阶的 Adams-Moulton 方法 Do 从步骤 i -> i+1 之前访问点，形式上是一个在 AB2
y[i+1] = y[i] + 0.5*dt * (3*f[i] - f[i-1])
对于第一个实现，以与 y 数组相同的方式声明 f 数组以逐字实现此公式是有意义的。仅保留一个 f 值的短数组可能更经济，这些值被移动或旋转以访问最后几个 f 值。
注意您需要使用其他类似或更高阶的方法初始化y[1]和f[1]。或者如果你想有一个 "pure" 运行 的方法，你需要初始化 y[-1] 和 f[-1] 并进一步返回以便 y[1] 可以用方法公式。