簡述

這里只考慮最為簡單的一種常微分方程

$$\frac{dy}{dx}=f(x,y)$$

然后這里的實例都是以下面這個方程來做展示的。

$$\frac{dy}{dx}=y?(1?y^2)$$
初值給定

$$y(0)=2$$

這個方程的精確解結(jié)果是下面這個方程

$$y(x)=\sqrt{\frac{4?e^{2x}}{4?e^{2x}?3}}$$

文章目錄

• • 簡述
  • 歐拉公式求解
  • • 簡單的理論推理
    • 代碼實現(xiàn)
    • • 實現(xiàn)后的效果
      • 代碼
      • 誤差畫圖
      • 誤差畫圖代碼
  • 改進(jìn)版歐拉公式
  • • 理解這個公式
    • 改進(jìn)版本的畫圖
    • 歐拉算法和改進(jìn)版歐拉算法的比較
    • • • 加上絕對值再來看
  • 累積誤差和分步的誤差
  • • 圖像
    • 代碼

歐拉公式求解

歐拉公式非常簡潔。(歐拉果然大佬！！！)

$$\begin{gathered} y_{n+1}=y_n+h?f(x_n,y_n) \\ {x}_n=x_0+n?h \end{gathered}$$

• h是步長

簡單的理論推理

其實非常直觀，將上面的第一個式子變形一下，有

$$\frac{dy}{dx}=f(x,y)=\frac{y_{n+1}?y_n}{h}$$

是不是非常熟悉，

微商的推導(dǎo)公式右邊的加上一個取極限（$h\to 0$）

之前的第一個式子，就是一階泰勒展開。

代碼實現(xiàn)

實現(xiàn)后的效果

x取值精確值數(shù)值逼近值誤差

0	2.0	2	0.0
0.1	1.6096571705090292	1.4	0.20965717050902932
0.2	1.4181045558702239	1.2656	0.15250455587022382
0.3	1.303667649645571	1.1894433603584	0.11422428928717099
0.4	1.2281238419433613	1.1401081630148027	0.08801567892855866
0.5	1.175177856120688	1.1059224047188059	0.06925545140188216
0.6	1.1365806251915203	1.0812532167953721	0.05532740839614814
0.7	1.1076620270914186	1.0629683037980915	0.0446937232933271
0.8	1.085560957285936	1.0491601738751934	0.03640078341074249
0.9	1.0684188538447474	1.0385912416407315	0.029827612204015974
1	1.0549729219451955	1.0304204608015664	0.02455246114362919

代碼

import numpy as np# dy/dx = f(x,y) # Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1 while x <= 1:loss = F(x) - yprint(x, '|', F(x), '|', y, '|', loss)y = y + h * f(y)x += 0.1

誤差畫圖

誤差畫圖代碼

import numpy as np import matplotlib.pyplot as plt# dy/dx = f(x,y) # Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1 losses = [] while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y = y + h * f(y)x += 0.1plt.plot(np.arange(0, 1.1, 0.1), losses) plt.plot(np.arange(0, 1.1, 0.1), losses, 'r*') plt.show()

改進(jìn)版歐拉公式

之前提到的歐拉公式簡單，容易實現(xiàn)，但是效率卻有點(diǎn)低。
所以有歐拉公式的改進(jìn)版本。

$$\begin{gathered} {\bar{y}}_{n+1}=y_n+h?f(x_n,y_n) \\ {y}_{n+1}=y_n+\frac{h}{2}?(f(x_n,y_n)+f(x_{n+1},{\bar{y}}_{n+1})) \end{gathered}$$

x的話，任然是按照步長移動的，所以，這里就不作贅述了。

理解這個公式

書上的話，是用定積分的方式推理出來的。
但是這個具體含義，也沒有給出什么東西來。
但是實際上，這個公式的邏輯含義是非常漂亮的！

將上面的兩個公式聯(lián)立起來。
$$\frac{y_{n+1}?y_n}{h}=\frac{1}{2}?(\frac{{\bar{y}}_{n+1}?y_n}{h}+f(x_{n+1},{\bar{y}}_{n+1}))$$

之前的歐拉方法是用這里的 $\frac{{\bar{y}}_{n+1}?y_n}{h}$來近似$\frac{dy}{dx}$
而現(xiàn)在我們逼近的其實是，用歐拉公式得到的導(dǎo)數(shù)，和實際上的在該點(diǎn)的導(dǎo)數(shù)值的均值。
從這個角度上來看，其實會是更貼近于正確的結(jié)果。

• 注意到： 我們這里用了$y_{n+1}$或者說是$y_n$是正確解。
• 用這個其實是合理的，因為，這里的算法是基于之前的歐拉算法的，而歐拉算法本來就是會近似逼近正確的結(jié)果。所以，$y_n$本身就是會逐漸接近于正確的結(jié)果。所以之前的假設(shè)是合理的。
• 而且，留心到，這是在原來的基礎(chǔ)上，做了新的加速的過程。

歐拉公式的預(yù)估值，和代入其的校正值的平均值，就是改進(jìn)歐拉方法的計算結(jié)果

改進(jìn)版本的畫圖

import numpy as np import matplotlib.pyplot as plt# dy/dx = f(x,y) # Advanced Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1 losses = [] while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1+y2)/2x += 0.1plt.plot(np.arange(0, 1.1, 0.1), losses) plt.plot(np.arange(0, 1.1, 0.1), losses, 'r*') plt.show()

歐拉算法和改進(jìn)版歐拉算法的比較

import numpy as np import matplotlib.pyplot as plt# dy/dx = f(x,y) # Advanced Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1def AdvancedEular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1 + y2) / 2x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y = y + h * f(y)x += 0.1return lossesadEl = AdvancedEular() El = Eular() plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Advance Eular') plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular') plt.plot(np.arange(0, 1.1, 0.1), El, 'b*') plt.legend() plt.show()

加上絕對值再來看

import numpy as np import matplotlib.pyplot as plt# dy/dx = f(x,y) # Advanced Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1def AdvancedEular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1 + y2) / 2x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y = y + h * f(y)x += 0.1return lossesadEl = AdvancedEular() El = Eular() plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Advance Eular') plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular') plt.plot(np.arange(0, 1.1, 0.1), El, 'b*') plt.legend() plt.show()

累積誤差和分步的誤差

前面的話，我們是用一個不那么正確的結(jié)果來放到公式中來推導(dǎo)出正確的結(jié)果。下面我們嘗試每一步都用正確的結(jié)果來推導(dǎo)出對應(yīng)y值。

x取值精確值數(shù)值逼近值(每步使用對應(yīng)的精確值推出)誤差

0	2.0	2	0.0
0.1	1.6096571705090292	1.4	0.20965717050902932
0.2	1.4181045558702239	1.35356132529304	0.06454323057718381
0.30000000000000004	1.303667649645571	1.274731273707409	0.028936375938162007
0.4	1.2281238419433613	1.2124696651611984	0.015654176782162965
0.5	1.175177856120688	1.1656997597866852	0.009478096334002872
0.6	1.1365806251915203	1.1303985272996449	0.006182097891875404
0.7	1.1076620270914186	1.103413438852542	0.004248588238876527
0.7999999999999999	1.085560957285936	1.082527495787655	0.003033461498280987
0.8999999999999999	1.0684188538447474	1.0661899261885104	0.0022289276562370564
0.9999999999999999	1.0549729219451955	1.0532987133870213	0.0016742085581742394

圖像

代碼

import numpy as np import matplotlib.pyplot as plt# dy/dx = f(x,y) # Euler formula f = lambda x: x * (1 - x ** 2) F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5 x = 0 y = 2 h = 0.1def stepEular():x = 0y = 2h = 0.1losses = []while x <= 1:ty = F(x)loss = ty - yprint(x, '|', ty, '|', y, '|', loss)losses.append(abs(loss))y = ty + h * f(ty)x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y = y + h * f(y)x += 0.1return lossesadEl = stepEular() El = Eular() plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Step Eular') plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular') plt.plot(np.arange(0, 1.1, 0.1), El, 'b*') plt.legend() plt.show()

總結(jié)

以上是生活随笔為你收集整理的常微分方程数值求解【python】的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。