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

1、精選優(yōu)質(zhì)文檔-傾情為你奉上解微分方程的歐拉法，龍格-庫(kù)塔法及其MATLAB簡(jiǎn)單實(shí)例歐拉方法（Euler method)用以對(duì)給定初值的常微分方程(即初值問(wèn)題)求解分為前進(jìn)EULER法、后退EULER法、改進(jìn)的EULER法。缺點(diǎn)：歐拉法簡(jiǎn)單地取切線的端點(diǎn)作為下一步的起點(diǎn)進(jìn)行計(jì)算，當(dāng)步數(shù)增多時(shí)，誤差會(huì)因積累而越來(lái)越大。因此歐拉格式一般不用于實(shí)際計(jì)算。改進(jìn)歐拉格式：為提高精度，需要在歐拉格式的基礎(chǔ)上進(jìn)行改進(jìn)。采用區(qū)間兩端的斜率的平均值作為直線方程的斜率。改進(jìn)歐拉法的精度為二階。算法為：微分方程的本質(zhì)特征是方程中含有導(dǎo)數(shù)項(xiàng)，數(shù)值解法的第一步就是設(shè)法消除其導(dǎo)數(shù)值。對(duì)于常微分方程： xa,by(

2、a) = y0可以將區(qū)間a,b分成n段，那么方程在第xi點(diǎn)有y'(xi) = f(xi,y(xi)，再用向前差商近似代替導(dǎo)數(shù)則為：在這里，h是步長(zhǎng)，即相鄰兩個(gè)結(jié)點(diǎn)間的距離。因此可以根據(jù)xi點(diǎn)和yi點(diǎn)的數(shù)值計(jì)算出yi+1來(lái)：i=0,1,2,L這就是向前歐拉格式。改進(jìn)的歐拉公式：將向前歐拉公式中的導(dǎo)數(shù)f(xi,yi)改為微元兩端導(dǎo)數(shù)的平均，即上式便是梯形的歐拉公式?？梢?jiàn)，上式是隱式格式，需要迭代求解。為了便于求解，使用改進(jìn)的歐拉公式： 數(shù)值分析中，龍格庫(kù)塔法（Runge-Kutta）是用于模擬常微分方程的解的重要的一類隱式或顯式迭代法。實(shí)際上，龍格-庫(kù)塔法是歐拉方法的一種推廣，向

3、前歐拉公式將導(dǎo)數(shù)項(xiàng)簡(jiǎn)單取為f(xn,yn),而改進(jìn)的歐拉公式將導(dǎo)數(shù)項(xiàng)取為兩端導(dǎo)數(shù)的平均。龍格-庫(kù)塔方法的基本思想：在區(qū)間xn,xn+1內(nèi)多取幾個(gè)點(diǎn)，將他們的斜率加權(quán)平均，作為導(dǎo)數(shù)的近似。 龍格庫(kù)塔法的家族中的一個(gè)成員如此常用，以至于經(jīng)常被稱為“RK4”或者就是“龍格庫(kù)塔法”。令初值問(wèn)題表述如下。則，對(duì)于該問(wèn)題的RK4由如下方程給出：其中這樣，下一個(gè)值(yn+1)由現(xiàn)在的值(yn)加上時(shí)間間隔(h)和一個(gè)估算的斜率的乘積決定。該斜率是以下斜率的加權(quán)平均：    k1是時(shí)間段開(kāi)始時(shí)的斜率；    k2是時(shí)間段中點(diǎn)的斜率，通過(guò)歐拉

4、法采用斜率k1來(lái)決定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)值：RK4法是四階方法，也就是說(shuō)每步的誤差是h5階，而總積累誤差為h4階。注意上述公式對(duì)于標(biāo)量或者向量函數(shù)(y可以是向量)都適用。  例子： 下面給出了數(shù)值求解該微分方程的簡(jiǎn)單程序。其中y1,y2,y3,y4分別為向前歐拉公式，改進(jìn)的歐拉公式，4級(jí)4階龍格-庫(kù)塔公式及精確解。h=0.1;x=0:h:1;y1=ze

5、ros(size(x);y1(1)=1;y2=zeros(size(x);y2(1)=1;y3=zeros(size(x);y3(1)=1;for i1=2:length(x)    y1(i1)=y1(i1-1)+h*(y1(i1-1)-2*x(i1-1)/y1(i1-1);    k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);    k2=y2(i1-1)+h*k1-2*x(i1)/(y2(i1-1)+h*k1);    y2(i1)=y2(i1-1)+h

6、*(k1+k2)/2;            k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);    k2=y2(i1-1)+h*k1/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k1/2);    k3=y2(i1-1)+h*k2/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k2/2);    k4=y2(i1-1)+h*k3-2*(x(i1-1)+h)/(y2(i1-1)+h*k3);    y3(i1)=y3(i1-1)+(k1+2*k2+2*k3+k4)*h/6;    endy4=sqrt(1+2*x);%plot(x,y1,x,y2,x,y3,x,y4)%legend('y1','y2','y3'