解微分方程歐拉法-R-K法及其MATLAB實(shí)例

1、歐拉方法(Euler method)用以對(duì)給定初值的常微分方程(即初值問題)求解 分為前進(jìn)EULER法、后退EULER法.改進(jìn)的EULER法。缺點(diǎn)：歐拉法簡單地取切線的端點(diǎn)作為下一步的起點(diǎn)進(jìn)行計(jì)算，當(dāng)步數(shù)增多時(shí)，誤差會(huì) 因積累而越來越大。因此歐拉格式一般不用于實(shí)際計(jì)算。改進(jìn)歐拉格式：為提髙精度，需要在歐拉格式的基礎(chǔ)上進(jìn)行改進(jìn)。釆用區(qū)間兩端的斜率的平均值 作為直線方程的斜率。改進(jìn)歐拉法的精度為二階。算法為：微分方程的本質(zhì)特征是方程中含有導(dǎo)數(shù)項(xiàng)，數(shù)值解法的第一步就是設(shè)法消除其導(dǎo) 數(shù)值。對(duì)于常微分方程：x丘a,by(a) = yO可以將區(qū)間a,b分成n段,那么方程在第xi點(diǎn)有y (xi) = f (

2、xi, y (xi), 再用向前差商近似代替導(dǎo)數(shù)則為：在這里，h是步長，即相鄰兩個(gè)結(jié)點(diǎn)間的距離。因此可以根據(jù)xi點(diǎn)和yi點(diǎn)的數(shù) 值計(jì)算出yi+1來：如=y-i + h x f (坯 yi)i=0,l,2,L這就是向前歐拉格式。改進(jìn)的歐拉公式：將向前歐拉公式中的導(dǎo)數(shù)f(xi,yi)改為微元兩端導(dǎo)數(shù)的平均，即上式便是梯形的歐拉公式?？梢姡鲜绞请[式格式，需要迭代求解。為了便于求解，使用改進(jìn)的歐拉公式: 數(shù)值分析中，龍格一庫塔法(Runge-Kutta)是用于模擬常微分方程的解的重要 的一類隱式或顯式迭代法。實(shí)際上，龍格-庫塔法是歐拉方法的一種推廣，向前歐拉公式將導(dǎo)數(shù)項(xiàng)簡單取為 f(xn,yn),而

3、改進(jìn)的歐拉公式將導(dǎo)數(shù)項(xiàng)取為兩端導(dǎo)數(shù)的平均。龍格-庫塔方法的基本思想：在區(qū)間xn,xn+l內(nèi)多取幾個(gè)點(diǎn)，將他們的斜率加權(quán)平均，作為導(dǎo)數(shù)的近似。龍格庫塔法的家族中的一個(gè)成員如此常用，以至于經(jīng)常被稱為RK4”或者就是 龍格庫塔法”。令初值問題表述如下。/ =y(io) = yo則，對(duì)于該問題的RK4由如下方程給出：h如+1 = yn + 6(局 + 2応2 + 2 應(yīng)3 + 魅)其中fcl = f (tn,yn)7 丄 hhT上2 = f (如+,羅幾+寸f 工 hhr上3 = f (尙+, +寸2也=f (tn + h,yn + hk3)這樣，下一個(gè)值()由現(xiàn)在的值W加上時(shí)間間隔(力)和一個(gè)估算的

4、斜率的乘積 決定。該斜率是以下斜率的加權(quán)平均：ki是時(shí)間段開始時(shí)的斜率;k2是時(shí)間段中點(diǎn)的斜率，通過歐拉法采用斜率kl來決定y在點(diǎn)tn + h/2的值；k3也是中點(diǎn)的斜率，但是這次采用斜率k2決定y值；k4是時(shí)間段終點(diǎn)的斜率，其y值用k3決定。當(dāng)四個(gè)斜率取平均時(shí)，中點(diǎn)的斜率有更大的權(quán)值：Ay + 2盪+ 2應(yīng)3 +魅 slope =.6RK4法是四階方法，也就是說每步的誤差是階，而總積累誤差為川階。注意上述公式對(duì)于標(biāo)量或者向量函數(shù)可以是向量)都適用。例子：下面給出了數(shù)值求解該微分方程的簡單程序。其中yI,y2,y3,y4分別為向前歐拉公式，改進(jìn)的歐拉公式，4級(jí)4階龍格-庫塔 公式及精確解。h二

5、；x=0：h：1；yl=zeros (size (x)；yl (1)=1;y2二zeros(size(x)；y2(l)=l;y3=zeros(size(x)；y3(l)=l;for i1=2：length(x) yl(il)=yl(il-l)+h*(yl(il-l)-2*x(il-l)/yl(il-l)； kl=y2(il-l)-2*x(il-l)/y2(il-l)；k2=y2(il-l)+h*kl-2*x(il)/(y2(il-l)+h*kl)； y2(il)=y2(il-l)+h* (kl+k2) /2；kl=y2(il-l)-2*x(il-l)/y2(il-l)；k2=y2(il-l)+h*k1/22*(x(ill)+h/2)/(y2(i1-1)+h*k1/2)； k3=y2(il-l)+h*k2/2-2*(x(il-1)+h/2)/(y2(il-l)+h*k2/2)； k4=y2(il-l) +h*k32*(x(il-l) +h)/(y2(il-l)+h*k3)； y3(il)=y3(il-l)+(kl+2*k2+2*k3+k4)*h/6；endy4=sqrt(l+2*x);%plot(x,yl,x,y2,x,y3,x,y4)