如何以相同的步长将欧拉方法与二阶龙格库塔方法进行比较？

Question

我有两种用于数值微分方程问题的算法，一种称为 Euler 方法，另一种称为二阶 Runge Kutta(RK2) 。本质上，欧拉方法和 RK2 近似于微分方程的解。唯一的区别是它们使用不同的公式（欧拉使用泰勒级数的一阶导数，而 RK2 是泰勒级数的二阶导数）。我正在尝试更正我编写的一些代码，以便它返回以下解决方案，但是，当涉及到 RK2 时，我的代码没有返回正确的值，而是返回以下内容，请注意，在我的解决方案中，我将步长 h 称为 dt。我在下面提供了用于创建此代码的代码，然后是二阶 Runge Kutta 方法的数值示例，该方法以数值方式工作。我有兴趣展示 RK2 的收敛速度比 Euler 的方法更快。

Answer 1

有两个小细节：

1. 在 RK-2 循环中，您使用的是 z，但要存储您使用的值 i

2. 在初始代码中，您在更新 K2 时使用了 y+K1/K2，这是错误的。我看到你在第二个代码中修复它。

所以，固定代码是：

import matplotlib.pyplot as plt

import numpy as np

from math import exp  # exponential function

dy = lambda x, y: x * y

f = lambda x: exp(x ** 2 / 2)  # analytical solution function

x_final = 2

# analytical solution

x_a = np.arange(0, x_final, 0.01)

y_a = np.zeros(len(x_a))

for i in range(len(x_a)):

    y_a[i] = f(x_a[i])

plt.plot(x_a, y_a, label="analytical")

# Container for step sizes dt /dt

dt = 0.5

x = 0

y = 1

print("dt = " + str(dt))

print("x \t\t y (Euler) \t y (analytical)")

print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_e = np.zeros(n + 1)

y_e = np.zeros(n + 1)

x_e[0] = x

y_e[0] = y

#Plot Euler's method

for i in range(n):

    y += dy(x, y) * dt

    x += dt

    print("%f \t %f \t %f" % (x, y, f(x)))

    x_e[i + 1] = x

    y_e[i + 1] = y

plt.plot(x_e, y_e, "x-", label="Euler dt=" + str(dt))

###################33

# Runge-Kutta's method 2'nd order (RK2)

x = 0

y = 1

print("dt = " + str(dt))

print("x \t\t y (rk2) \t y (analytical)")

print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_r = np.zeros(n + 1)

y_r = np.zeros(n + 1)

x_r[0] = x

y_r[0] = y

# Plot the RK2

for i in range(n):

    K1 = dt*dy(x,y) # Step 1

    K2 = dt*dy(x+dt/2,y+K1/2) # Step 2

    y += K2 # Step 3

    x += dt

    print("%f \t %f \t %f" % (x, y, f(x)))

    x_r[i + 1] = x

    y_r[i + 1] = y

plt.plot(x_r, y_r, "s-", label="RK2 dt=" + str(dt))

plt.title("numerical differential equation problem")

plt.xlabel("x")

plt.ylabel("y")

plt.legend()

plt.show()