我有兩種用于數(shù)值微分方程問題的算法，一種稱為 Euler 方法，另一種稱為二階 Runge Kutta(RK2) 。本質(zhì)上，歐拉方法和 RK2 近似于微分方程的解。唯一的區(qū)別是它們使用不同的公式（歐拉使用泰勒級數(shù)的一階導(dǎo)數(shù)，而 RK2 是泰勒級數(shù)的二階導(dǎo)數(shù)）。我正在嘗試更正我編寫的一些代碼，以便它返回以下解決方案，但是，當(dāng)涉及到 RK2 時(shí)，我的代碼沒有返回正確的值，而是返回以下內(nèi)容，請注意，在我的解決方案中，我將步長 h 稱為 dt。我在下面提供了用于創(chuàng)建此代碼的代碼，然后是二階 Runge Kutta 方法的數(shù)值示例，該方法以數(shù)值方式工作。我有興趣展示 RK2 的收斂速度比 Euler 的方法更快。
1 回答
青春有我
TA貢獻(xiàn)1784條經(jīng)驗(yàn) 獲得超8個(gè)贊
有兩個(gè)小細(xì)節(jié)：
在 RK-2 循環(huán)中，您使用的是 z，但要存儲您使用的值 i
在初始代碼中，您在更新 K2 時(shí)使用了 y+K1/K2，這是錯誤的。我看到你在第二個(gè)代碼中修復(fù)它。
所以，固定代碼是：
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()
添加回答
舉報(bào)
0/150
提交
取消