解微分方程歐拉法RK法及其MATLAB實例

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

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

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

4、tn + h/2的值；    k3也是中點的斜率，但是這次采用斜率k2決定y值；    k4是時間段終點的斜率，其y值用k3決定。當(dāng)四個斜率取平均時，中點的斜率有更大的權(quán)值：RK4法是四階方法，也就是說每步的誤差是h5階，而總積累誤差為h4階。注意上述公式對于標(biāo)量或者向量函數(shù)(y可以是向量)都適用。  例子： 下面給出了數(shù)值求解該微分方程的簡單程序。其中y1,y2,y3,y4分別為向前歐拉公式，改進(jìn)的歐拉公式，4級4階龍格-庫塔公式及精確解。h=0.1;x=0:h:1;y1=zeros(size(x);y

5、1(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*(k1+k2)/2;&#

6、160;           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','