欧拉法逼近太准确

Question

我正在尝试编写欧拉逼近法的代码来解决谐振子问题，然后在绘图上将其与精确解进行比较。然而，我发现我的近似值似乎太好了（即使对于更大的步长），尤其是因为我知道误差应该随着时间的推移而增加。我在我的欧拉方法代码中做错了什么吗？都写在下面了。

import numpy as np
import matplotlib.pyplot as plt

w = 1.0 #frequency of oscillation
x_0 = 1.0 #initial position
h = 0.6 #step size
t_0 = 0 #starting time
n = 20
t = np.arange(t_0, n, h) #t = t_0 + nh (h has units of seconds)
y = np.array([[0.0, 0.0]]*len(t))

def x(t, x_0, w): #exact solution
    return x_0*np.cos(w*t)
def x_prime(t, x_0, w): #first derivative of solution (which is also v)
    return -x_0*w*np.sin(w*t)

y[0] = [x_0, 0.0] #starting velocity is 0
for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + np.array([x_prime(t[i-1], x_0, w), -w**2 * x(t[i-1], x_0, w)]) * h

xlist = [item[0] for item in y] #list of x values at each t value
vlist = [item[1] for item in y] #list of velocity values at each t value

plt.plot(t, xlist, label='Euler approximation')
plt.plot(t, x(t, x_0, w), label='Exact solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend(loc='upper right')

我的这段代码输出为：

Answer 1

您的问题是您根据精确解计算欧拉方法的斜率。您应该改用之前的近似值，

# Euler's method
y[i] = y[i-1] + np.array([y[i-1,1], -w**2 * y[i-1,0]]) * h

或使用函数

def f(t,y): x, dotx = y; return np.array([ dotx, -w**2*x])

...

for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + f(t[i-1],y[i-1]) * h

然后您将获得预期的快速发散行为，此处为 h=0.1：

顺便说一句，写起来可能更简单

xlist, vlist = y.T

欧拉法逼近太准确

Euler method approximation is too accurate

python

python-2.7

numerical-methods

approximation