龙格-库塔四次方法

Question

我正在实现Runge–Kutta四阶方法来求解两个方程组。在

h是分段数，所以T/h是步长。在

def cauchy(f1, f2, x10, x20, T, h):
    x1 = [x10]
    x2 = [x20]

    for i in range(1, h):
        k11 = f1((i-1)*T/h, x1[-1], x2[-1])
        k12 = f2((i-1)*T/h, x1[-1], x2[-1])
        k21 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12)
        k22 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12)
        k31 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22)
        k32 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22)
        k41 = f1((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32)
        k42 = f2((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32)

        x1.append(x1[-1] + T/h/6*(k11 + 2*k21 + 2*k31 + k41))
        x2.append(x2[-1] + T/h/6*(k12 + 2*k22 + 2*k32 + k42))

    return x1, x2

然后我在这个系统上测试它：

^{pr2}$

它似乎运行良好（我还用不同的初始值和不同的函数对其进行了测试：一切正常）：

x10 = 1
x20 = 1
T = 1
h = 10

x1, x2 = cauchy(f1, f2, x10, x20, T, h)
t = np.linspace(0, T, h)

plt.xlabel('t')
plt.ylabel('x1')
plt.plot(t, true_x1(t), "blue", label="true_x1")
plt.plot(t, x1, "red", label="approximation_x1")
plt.legend(bbox_to_anchor=(0.97, 0.27))
plt.show()

plt.xlabel('t')
plt.ylabel('x2')
plt.plot(t, true_x2(t), "blue", label="true_x2")
plt.plot(t, x2, "red", label="approximation_x2")
plt.legend(bbox_to_anchor=(0.97, 0.97))
plt.show()

然后我想检查错误是否是O(step^4)的顺序，所以我减少了步骤和计算错误，如下所示：

step = []
x1_error = []
x2_error = []
for segm in reversed(range(10, 1000)):
    x1, x2 = cauchy(f1, f2, x10, x20, T, segm)
    t = np.linspace(0, T, segm)
    step.append(1/segm)
    x1_error.append(np.linalg.norm(x1 - true_x1(t), np.inf))
    x2_error.append(np.linalg.norm(x2 - true_x2(t), np.inf))

我明白了：

plt.plot(step, x1_error, label="x1_error")
plt.plot(step, x2_error, label="x2_error")
plt.legend()

所以，误差是线性的。这真的很奇怪，因为它应该是O(step^4)的顺序。谁能告诉我我做错了什么吗？在

Answer 1

for i in range(1, h):

这将从1迭代到h-1。由于缺少最后一步，从时间T-T/h到时间T的精确解的差是O(T/h)。在

因此使用

^{pr2}$

对于从i-1到{}的h步骤，或

for i in range(h):

对于h步骤，从i到{}。在

---

另外，np.linspace(0,1,4)将产生4等间距数字，其中第一个是0，最后一个是{}，结果是

array([ 0.        ,  0.33333333,  0.66666667,  1.        ])

这可能不是你所期望的。所以用上面的修正使用

t = np.linspace(0, T, segm+1)

在两个计算中使用相同的时间点。在

---

如果使用通常意义上的字母，那么遵循代码会更容易，其中h或{}是步长，N是步骤数。然后在循环h=T/N或dt=T/N之前定义，以避免在函数调用中重复使用T/N。在