數(shù)值分析2常微分方程解

常微分方程的數(shù)值解法電子科技大學(xué)常微分方程的數(shù)值解引言簡單的數(shù)值方法Runge-Kutta方法一階常微分方程組和高階方程在高等數(shù)學(xué)中我們見過以下常微分方程：6.1引言(1)，(2)式稱為初值問題，(3)式稱為邊值問題。在實(shí)際應(yīng)用中還經(jīng)常需要求解常微分方程組：本章主要研究問題(1)的數(shù)值解法，對(2)~(4)只作簡單介紹。（其中L為Lipschitz常數(shù)）則初值問題（1）存在唯一的連續(xù)解。

考慮一階常微分方程初值問題其中，y=y(x)是未知函數(shù)，y(x0)=y0

是初值條件，而f(x,y)是給定的二元函數(shù).由常微分方程理論知，若f(x)在x[a,b]連續(xù)且f滿足對y的Lipschitz條件：常微分方程的數(shù)值解法有單步法和多步法之分：單步法：在計(jì)算yn＋1時(shí)只用到前一點(diǎn)yn

的值；多步法：計(jì)算yn＋1時(shí)不僅利用yn，還要利用yn-1,yn-2,…….，一般k步法要用到y(tǒng)n,yn-1,yn-2,…….,yn-k+1。求問題（1）的數(shù)值解，就是要尋找解函數(shù)在一系列離散節(jié)點(diǎn)x1<x2<……<xn<xn+1

上的近似值y1,y

2,…,yn

。為了計(jì)算方便，可取xn=x0+nh，(n=0,1,2,…)，h稱為步長。6.2簡單的數(shù)值方法一、歐拉（Euler）方法在x=x0處，用差商代替導(dǎo)數(shù)：由得同理，在x=xn

處，用差商代替導(dǎo)數(shù)：由得若記則上式可記為此即為求解初值問題的Euler方法，又稱顯式Euler方法。Pn+1yOxx0x1x2xnP0P1P2Pny=f(x)xn+1Euler方法的幾何意義：（Euler折線法）例:用Euler方法求解常微分方程初值問題并將數(shù)值解和該問題的解析解比較。解：Euler方法的具體格式：

xn

y(xn) yn

yn-y(xn) 0.0 0 0 00.2 0.1923 0.2000 0.00770.4 0.3448 0.3840 0.03920.6 0.4412 0.5170 0.07580.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.09241.2 0.4918 0.5705 0.07871.4 0.4730 0.5354 0.0624取h=0.2,xn=nh,(n=0,1,2…,15),f(x,y)=y/x–2y2

計(jì)算中取f(0,0)=1.計(jì)算結(jié)果如下：xn

y(xn) yn

yn-y(xn)1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.03592.0 0.4000 0.4268 0.02682.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.01472.6 0.3351 0.3459 0.01082.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057由表中數(shù)據(jù)可以看到，微分方程初值問題的數(shù)值解和解析解的誤差一般在小數(shù)點(diǎn)后第二位或第三位小數(shù)上，這說明Euler方法的精度是比較差的。O:數(shù)值解—:準(zhǔn)確解h=0.2;y(1)=0.2;x=h:h:3;forn=1:14xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0:h:3;y0=x0./(1+x0.^2);plot(x0,y0,x,y,x,y,'o')若直接對y`=f(x,y)在[xn,xn+1]積分，利用數(shù)值積分中的左矩形公式：此即為Euler公式。設(shè)y(xn)=yn，則得若用右矩形公式：得上式稱后退的Euler方法，又稱隱式Euler方法?？捎玫ㄇ蠼猓撼踔担旱簁=0,1,……因故當(dāng)hL<1時(shí)，迭代法收斂。二、梯形方法由利用梯形求積公式：得上式稱梯形方法，是一種隱式方法。用迭代法求解：初值：迭代：k=0,1,……因故當(dāng)hL/2<1時(shí)，迭代法收斂。由以上分析可以看出，隱式方法的計(jì)算比顯式方法復(fù)雜，需要用迭代法求解非線性方程才能得出計(jì)算結(jié)果?？刹捎脤@式Euler格式與梯形格式結(jié)合使用的方法來避免求解非線性方程。記再用梯形格式計(jì)算：－－預(yù)測－－校正上面兩式統(tǒng)稱預(yù)測－校正法，又稱改進(jìn)的Euler方法。三、單步法的局部截?cái)嗾`差和精度單步法的一般形式為：（

與f有關(guān)）顯式單步法形式為：整體截?cái)嗾`差：從x0開始，考慮每一步產(chǎn)生的誤差，直到xn，則有誤差稱為數(shù)值方法在節(jié)點(diǎn)xn處的整體截?cái)嗾`差。但en不易分析和計(jì)算，故只考慮從xn到xn+1的局部情況。定義：設(shè)y(x)是初值問題(1)的精確解，則稱為顯式單步法在節(jié)點(diǎn)xn+1處的局部截?cái)嗾`差。若存在最大整數(shù)p使局部截?cái)嗾`差滿足則稱顯式單步法具有p階精度或稱p階方法。注：將Tn+1表達(dá)式各項(xiàng)在xn處作Taylor展開，可得具體表達(dá)式。Euler方法的局部截?cái)嗾`差：故Tn+1=O(h2)，p=1，（設(shè)yn=y(xn)）其中稱局部截?cái)嗾`差主項(xiàng)。即Euler方法具1階精度。（設(shè)yn=y(xn)）故Tn+1=O(h3)，p=2，梯形方法的局部截?cái)嗾`差：局部截?cái)嗾`差主項(xiàng)為：梯形方法具2階精度。6.3Runge－Kutta方法一、Runge-Kutta方法的基本思想由Taylor展式Tn+1=O(hp+1)，若提高p，可提高精度。但因……高階導(dǎo)數(shù)計(jì)算復(fù)雜，故可從另外角度考慮。分析Euler公式及改進(jìn)的Euler公式：局部截?cái)嗾`差：O(h2)局部截?cái)嗾`差：O(h3)可用f(x,y)在某些點(diǎn)處值的線性組合得yn+1，增加計(jì)算f(x,y)的次數(shù)可提高階數(shù)。設(shè)法計(jì)算f(x,y)在某些點(diǎn)上的函數(shù)值，然后對這些函數(shù)值作線性組合，構(gòu)造近似計(jì)算公式，再把近似公式和解的泰勒展開式相比較，使前面的若干項(xiàng)吻合，從而獲得達(dá)到一定精度的數(shù)值計(jì)算公式。Runge-Kutta方法的基本思想：設(shè)ci,

λi,

μij為待定常數(shù)。上面第一個(gè)式子的右端在(xn,yn)作泰勒展開后,按h的冪次作升序排列：再與初值問題的精確解y(x)在點(diǎn)x=xn處的泰勒展開式相比較，使其有盡可能多的項(xiàng)重合。例如，要求就得到p個(gè)方程,從而定出參數(shù)ci，

i，

ij，再代入K1，K2，…，Kr的表達(dá)式，就可得到計(jì)算微分方程初值問題的數(shù)值計(jì)算公式：若

Tn+1=O(hp+1)，則稱其為p階r級(jí)R－K方法。上式稱為r級(jí)Runge－Kutta方法的計(jì)算公式。當(dāng)r=1時(shí)，就是Euler方法。要使Runge-Kutta公式具有更高的階p，就要增加r的值。下面我們只就r=2推導(dǎo)R-K方法。二、二階Runge-Kutta方法其中c1,c2,

2,

21待定。上式的局部截?cái)嗾`差為：又由利用二元函數(shù)的Taylor展開，得代入Tn+1的表達(dá)式，得即c1=1-a，

2=

21=1/(2a)要使上式p=2階，則需方程組解不唯一，可令c2=a

0

，則

滿足上述條件的公式都為2階R-K公式。稱中點(diǎn)公式，相當(dāng)于數(shù)值積分的中矩形公式：如取a=?，則c1=

c2=?，

2=

21=1，即為改進(jìn)Euler公式。若取a=1，則c1=0，c2=1，

2=

21

=?，得例:蛇形曲線的初值問題令f(x,y)=y/x–2y2,取f(0,0)=1,h=0.2,xn=hn,(n=1,2,…,15)2階龍格-庫塔公式計(jì)算格式：

k1=yn/xn–2yn

2,k2=(yn+hk1)/(xn+h)–2(yn+hk1)2yn+1=yn+0.5h[k1+k2]x0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2;y(1)=y0+.5*h*(k1+k2);forn=1:14k1=y(n)/x(n)-2*y(n)^2;k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2;y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.^2);subplot(221)plot(x,y1,x,y1,'b*')subplot(222)plot(x,y,x,y,'o')subplot(223)plot(x,y,x,y,'o',x,y1,x,y1,'b*')

三、三階與四階Runge-Kutta方法當(dāng)r=3時(shí)，R-K公式表示為其中為8個(gè)待定常數(shù)。上式的局部截?cái)嗾`差為類似二階的推導(dǎo)過程，將K2,K3按二元函數(shù)展開，使Tn+1=O(h4)，得方程有8個(gè)未知數(shù)，解不唯一。滿足該條件的公式統(tǒng)稱為三階R-K公式。其中一個(gè)常用公式為：當(dāng)r=4時(shí)，利用相同的推導(dǎo)過程，經(jīng)過較復(fù)雜的計(jì)算，可以得出四階R-K公式的成立條件。下列經(jīng)典公式是其中常用的一個(gè)：k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2),k4=f(xn+h,yn+hk3)

yn+1=yn+h[k1+2k2+2k3+k4]/6(n=0,1,2,·····,N)四階龍格-庫塔公式數(shù)值求解:狐貍和野兔問題——常微分方程組在一個(gè)生物圈中,只有野兔和狐貍兩種動(dòng)物,記y1

為野兔數(shù)量,y2

為狐貍數(shù)量.設(shè)有足夠的食物供野兔生存,而狐貍只以野兔為食物.

模型如下自變量x∈(0,15),初值條件為:y1(0)=20,y2(0)=20t0=0;t1=15; %設(shè)置區(qū)域y0=[2020]; %給定初值[t,y]=ode23(‘lot1’,t0,t1,y0); %求解plot(t,y) %繪圖functionz=lot1(x,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);求解常微分方程:(1)定義函數(shù);(2)調(diào)用ode236.5一階常微分方程組和高階微分方程的數(shù)值解法簡介一、一階常微分方程組的數(shù)值解法：下