如何将欧拉方法与步长相同的二阶龙格库塔方法进行比较?
How to compare Euler's method to a second-Order Runge-Kutta Method at the same stepsize?
对于数值微分方程问题,我有两种算法,一种称为欧拉方法,另一种称为二阶龙格库塔(RK2)。本质上,欧拉方法和 RK2 逼近微分方程的解。唯一的区别是它们使用不同的公式(欧拉使用泰勒级数的一阶导数,而 RK2 是泰勒级数的二阶导数)。
我正在尝试更正我编写的一些代码,以便 return 以下解决方案,
但是,当涉及到 RK2 时,我的代码没有 return 正确的值,而是 return 以下,
请注意,在我的解决方案中,我将步长 h 称为 dt。我在下面提供了我用来创建它的代码,然后是一个以数字方式工作的二阶 Runge Kutta 方法的数字示例。我有兴趣证明 RK2 的收敛速度比欧拉的方法快。
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 步放在允许我们创建 RK2 的 Euler 方法的更改上。
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 将其绘制在图表上
有两个小细节:
- 在 RK-2 循环中您使用 z,但要存储您使用 i 的值
- 在初始代码中,您在更新 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()
对于数值微分方程问题,我有两种算法,一种称为欧拉方法,另一种称为二阶龙格库塔(RK2)。本质上,欧拉方法和 RK2 逼近微分方程的解。唯一的区别是它们使用不同的公式(欧拉使用泰勒级数的一阶导数,而 RK2 是泰勒级数的二阶导数)。
我正在尝试更正我编写的一些代码,以便 return 以下解决方案,
但是,当涉及到 RK2 时,我的代码没有 return 正确的值,而是 return 以下,
请注意,在我的解决方案中,我将步长 h 称为 dt。我在下面提供了我用来创建它的代码,然后是一个以数字方式工作的二阶 Runge Kutta 方法的数字示例。我有兴趣证明 RK2 的收敛速度比欧拉的方法快。
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 步放在允许我们创建 RK2 的 Euler 方法的更改上。
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 将其绘制在图表上有两个小细节:
- 在 RK-2 循环中您使用 z,但要存储您使用 i 的值
- 在初始代码中,您在更新 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()