一个多元微分方程的数值求解

Question

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

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

初始条件为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

我认为问题在于你对欧拉方案的理解。这是一个非常简单的常微分方程，事实上它是谐波振荡系统的教科书上的例子。欧拉格式基本上基于N=T/（dt），其中N是步数，T是最终时间，dt是步长。因此，如果步长或最终时间相对于N不小，则解会向上漂移

不需要使用RK4，Euler将完成工作。但诀窍是，你需要一个小dt。我已经以更清晰的方式重写了你的方案，并分别使用了适当的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()