python - 使用 Numpy 的中心差异算法

Question

我正在尝试使用 Numpy 实现一个相当简单的算法来数值求解微分方程组。我使用所谓的欧拉前向方法很好地工作（即解决方案是预期的指数衰减），但它使用中心差异方案产生有趣的结果（参见下图中的红色和绿色曲线，与橙色和蓝色相比） .

我怀疑这可能与 Numpy（以及更普遍的 Python）处理分配的方式有关（是否只是一个名称，是否通过引用过去，是否np.copy()应该使用切片表示法或 a 等）。老实说，我已经阅读了几篇关于这个主题的文章，但它仍然让我感到困惑（我是 Fortran 的新手）。

我正在使用并重现上图的代码的简化版本如下。

import numpy as np
import matplotlib.pyplot as plt


N = 100
Dt = 10 / N
t = np.linspace(0, 10, N+1)

y1 = np.zeros(N+1)
y2 = np.zeros(N+1)

##########################################
# problem definition (the result is an exponential decay)
y1[0] = 1
y2[0] = 2

def f(y):
    return -0.1 * y
##########################################



fig, ax = plt.subplots()

##########################################
# Euler-forward
y = np.ones(2)
y[0] = y1[0]
y[1] = y2[0]

for i in range(1, N+1):
    y = y + f(y) * Dt
    y1[i] = y[0]
    y2[i] = y[1]

ax.plot(t, y1)
ax.plot(t, y2)
##########################################


##########################################
# central differences
y = np.ones(2)
y[0] = y1[0]
y[1] = y2[0]

# first step is Euler-forward
y_new = y + f(y) * Dt
y_old = y
y = y_new

y1[1] = y[0]
y2[1] = y[1]

for i in range(2, N):
    f = np.cos(y)
    y_new = y_old + 2 * f * Dt
    y_old = y
    y = y_new
    y1[i] = y[0]
    y2[i] = y[1]

ax.plot(t, y1)
ax.plot(t, y2)
########################################## 

标签： pythonalgorithmnumpy