如何在相同步长下比较Euler方法和二阶RungeKutta方法？

Question

对于数值微分方程问题，我有两种算法，一种叫做Euler方法，另一种叫做二阶Runge-Kutta（RK2）。本质上，欧拉方法和RK2近似于微分方程的解。唯一的区别是它们使用不同的公式（欧拉使用泰勒级数的一阶导数，而RK2是泰勒级数的二阶导数）

我试图更正我编写的一些代码，以便它返回以下解决方案，

但是，我的代码在涉及RK2时没有返回正确的值，而是返回以下值：

注意，在我的解决方案中，我将步长h称为dt。我在下面提供了我用来创建它的代码，后面是一个二阶龙格-库塔方法的数值示例，该方法在数值上工作。我感兴趣的是证明RK2比Euler方法的收敛速度更快

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_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[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_n[i + 1] = x
    y_n[i + 1] = y

plt.plot(x_n, y_n, "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_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y

# Plot the RK2
for z in range(n):
    K1 = dt*dy(x,y) # Step 1
    K2 = dt*dy(x+dt/2,y+K1/K2) # Step 2
    y += K2 # Step 3
    x += dt
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_n[i + 1] = x
    y_n[i + 1] = y

plt.plot(x_n, y_n, "x-", label="RK2 dt=" + str(dt))

plt.title("Euler's method ")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

这是适用于RK2的数字代码。我已经将步骤1、2和3放在了Euler方法的更改上，该方法允许我们创建RK2

from math import exp
dy = lambda x,y: x*y
f = lambda x: exp(x**2/2)
x = 0
xn = 2
y = 1
dt = 0.5
n = int((xn)/dt)
print ('x \t\t y (RK2) \t y (analytical)')
print ('%f \t %f \t %f'% (x,y,f(x)))
# main loop
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)))

对不起，如果我问得不好。我是python新手，所以我仍在努力理解如何着手解决这些类型的问题。我的问题的总结是如何修复上述函数，以便计算正确的估计值，然后使用python中的matplotlib将其绘制在图形上

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()