亚洲在线久爱草,狠狠天天香蕉网,天天搞日日干久草,伊人亚洲日本欧美

為了賬號安全,請及時綁定郵箱和手機立即綁定
已解決430363個問題,去搜搜看,總會有你想問的

如何以相同的步長將歐拉方法與二階龍格庫塔方法進行比較?

如何以相同的步長將歐拉方法與二階龍格庫塔方法進行比較?

元芳怎么了 2022-10-25 15:15:23
我有兩種用于數值微分方程問題的算法,一種稱為 Euler 方法,另一種稱為二階 Runge Kutta(RK2) 。本質上,歐拉方法和 RK2 近似于微分方程的解。唯一的區別是它們使用不同的公式(歐拉使用泰勒級數的一階導數,而 RK2 是泰勒級數的二階導數)。我正在嘗試更正我編寫的一些代碼,以便它返回以下解決方案,但是,當涉及到 RK2 時,我的代碼沒有返回正確的值,而是返回以下內容,請注意,在我的解決方案中,我將步長 h 稱為 dt。我在下面提供了用于創建此代碼的代碼,然后是二階 Runge Kutta 方法的數值示例,該方法以數值方式工作。我有興趣展示 RK2 的收斂速度比 Euler 的方法更快。
查看完整描述

1 回答

?
青春有我

TA貢獻1784條經驗 獲得超8個贊

有兩個小細節:

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


查看完整回答
反對 回復 2022-10-25
  • 1 回答
  • 0 關注
  • 134 瀏覽
慕課專欄
更多

添加回答

舉報

0/150
提交
取消
微信客服

購課補貼
聯系客服咨詢優惠詳情

幫助反饋 APP下載

慕課網APP
您的移動學習伙伴

公眾號

掃描二維碼
關注慕課網微信公眾號