在没有 odeint 的情况下数值求解多元微分方程

Question

我很好奇。我知道这可以通过使用 odeint 来解决，但我正在尝试从头开始，我遇到了一个有趣的行为。

假设一个简单的振荡器，方程式 m * x_ddot + k * x = 0。 Wikipedia

初始条件为x0 != 0；

理论上，解是一个正弦函数。

但是，在 python 上，解是振幅不断增长的正弦波。现在我很好奇为什么会这样，因为它不应该发生。它与数值稳定性或类似的东西有关吗？比如，为什么会出现分歧？从物理学的角度来看，没有理由应该如此，那么为什么它会这样呢？

这是代码。

dt = 0.05
t_start = 0
t_finish = 20
t = 0
x1 = 1
X1 = []

x2 = 0
X2 = []

while t <= t_finish:    
    
    X1.append(x1)
    X2.append(x2)
    
    # state space representation
    x1_dot = x2
    x2_dot = -9.81*x1
    
    x1 += dt*x1_dot
    x2 += dt*x2_dot
    
    t += dt

# to make sure the vectors are of equal size for plotting
if len(X1) > len(time):
    X1 = X1[:len(X1)-1]
elif len(X1) < len(time):
    time = time[:len(time)-1]
plt.figure(figsize=(10,10))
plt.plot(time,X1)
plt.grid()

这是情节。

非常感谢你们提供的任何见解。

Answer 1

我觉得问题出在你对欧拉方案的理解上。这是一个非常简单的 ODE，实际上它是谐波振荡系统的教科书示例。欧拉方案几乎基于 N = T/(dt)，其中 N 是步数，T 是最终时间，dt 是步长。所以，如果你的步长或最终时间相对于 N 不小，解就会向上漂移。

不用RK4，Euler搞定。不过诀窍是你需要一个小的 dt。我把你的scheme改写得更清楚了，分别用了合适的dt。

import numpy as np
import matplotlib.pyplot as plt

# Parameters
t_finish = 20.0
dt = 0.00005 # Very small dt (infinitesmal dt)
n = int(t_finish/dt)

tvalue = np.linspace(0, t_finish, n+1)
x1 = np.zeros(n+1)
x2 = np.zeros(n+1)

# create canvas
plt.figure(figsize=(10, 5))

# Initialize
x1[0] = 1.0 # Not at zero
x2[0] = 0.0

# Simulation with Euler scheme
for i in range(n):
        t = (i+1)*dt
        x1[i+1] = x1[i] + x2[i]*dt
        x2[i+1] = x2[i] -9.81*x1[i]*dt
        
# Plot paths
plt.plot(tvalue, x1, label=r'$x_1$')
plt.plot(tvalue, x2, label=r'$x_2$')

# Add legend and axes labels
plt.legend(loc=0)
plt.xlabel(r'$t$')
plt.ylabel(r'$x_{t}$')
plt.show()

在没有 odeint 的情况下数值求解多元微分方程

Numerically Solving A Multivariate Differential Equation Without odeint

python

physics

numeric

numerical-methods

differential-equations