1. 程式人生 > >應用數學課堂筆記——常微分方程數值解

應用數學課堂筆記——常微分方程數值解

此為應用數學第5次課。

常微分方程數值方法

差分法、有限元法、譜方法等。這裡只介紹顯式尤拉法。

尤拉法及其變種

問題描述:在x[a,b]x \in [a,b]求解yy,滿足 y=f(x,y),y(x0)=y0 y'= f(x,y), y(x_0)=y_0

[a,b][a,b]等分成NN份,每份長度為hhx0=a,xi=a+ihx_0=a, x_i=a+ih。令y(xi)y(x_i)xix_i處的真實值,yiy_i為擬合值。

前向尤拉法

由泰勒展開,有 y(xi

)=y(xi+Δx)y(xi)Δx+o(Δx) y'(x_i) = \frac{y( x_i+\Delta x ) - y(x_i)}{\Delta x} + o(\Delta x)

近似得到, y(xi)y(xi+h)y(xi)h y'(x_i) \approx \frac{y(x_i+h)-y(x_i)}{h}

因為y(xi)=f(xi,y(xi))y'(x_i) = f(x_i, y(x_i))

)=f(xi,y(xi)),所以有 y(xi+1)y(xi)hf(xi,y(xi)) y(x_{i+1}) - y(x_i) \approx h f(x_i, y(x_i))

從而得到前向尤拉法, yi+1=yi+hf(xi,yi) y_{i+1} = y_i + h f(x_i, y_i)

後向尤拉法

類似的,有 y(xi+1)y(xi+1)y(xi)h y'(x_{i+1}) \approx \frac{y(x_{i+1})-y(x_i)}{h}

xi+1)hy(xi+1)y(xi)

由於y(xi+1)=f(xi+1,yi+1)y'(x_{i+1}) = f(x_{i+1}, y_{i+1}),得到 yi+1=yi+hf(xi+1,yi+1) y_{i+1} = y_i + h f(x_{i+1}, y_{i+1})

但由於yi+1y_{i+1}不預先知道,所以需要方程解出yi+1y_{i+1},為隱式尤拉法了。可以使用預估校正法來利用後向尤拉的思想。

預估校正法

先用前向尤拉法得到對yi+1y_{i+1}的估計yˉi+1\bar{y}_{i+1},即 yˉi+1=yi+hf(xi,yi) \bar{y}_{i+1} = y_i + h f(x_i, y_i)

再用後向尤拉法進行校正 yi+1=yi+hf(xi+1,yˉi+1) y_{i+1} = y_i + h f(x_{i+1}, \bar{y}_{i+1})

梯形法

綜合一下前向尤拉和後向尤拉,可以得到一個平均主義的變種: yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1)] y_{i+1} = y_i + \frac{h}{2} [ f(x_i, y_i) + f(x_{i+1}, y_{i+1}) ]

當然這仍然是一種隱式尤拉法。

在預估校正法中也可以採用類似的思想,有 yi+1=yi+h2[f(xi,yi)+f(xi+1,yˉi+1)] y_{i+1} = y_i + \frac{h}{2} [ f(x_i, y_i) + f(x_{i+1}, \bar{y}_{i+1}) ]

2階龍格庫塔

2階龍格庫塔如下: yi+1=yi+h[12k1+12k2]k1=f(xi,yi)k2=f(xi,yi+hk1) \begin{aligned} y_{i+1} &= y_i + h[\frac{1}{2} k_1 + \frac{1}{2} k_2] \\ k_1 &= f(x_i, y_i) \\ k_2 &= f(x_i, y_i + h k_1) \end{aligned}

可以看到,2階龍格庫塔和預估校正的梯形法非常類似,它們都使用了yˉi+1\bar{y}_{i+1}處的梯度作為