階RungeKutta法求解一階常微分方程

1、MATLAB語言及應(yīng)用大作業(yè)姓名：學(xué)號：學(xué)院：班級：題目編號：2013年10月134階Runge-Kutta法求解一階常微分方程。一、 Runge-Kutta法的數(shù)學(xué)理論龍格-庫塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進(jìn)行抑制，所以其實(shí)現(xiàn)原理也較復(fù)雜。該算法是構(gòu)建在數(shù)學(xué)支持的基礎(chǔ)之上的。龍格庫塔方法的理論基礎(chǔ)來源于泰勒公式和使用斜率近似表達(dá)微分，它在積分區(qū)間多預(yù)計算出幾個點(diǎn)的斜率，然后進(jìn)行加權(quán)平均，用做下一點(diǎn)的依據(jù)，從而構(gòu)造出了精度更高的數(shù)值積分計算方法。如果預(yù)先求兩個點(diǎn)的斜率就是二階龍格庫塔法，如果預(yù)先取四個點(diǎn)就是四階龍格庫塔法

2、。一階常微分方程可以寫作：y'=f(x,y),使用差分概念。(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，極限為Yn'）Yn+1=Yn+h*f(Xn,Yn)另外根據(jù)微分中值定理，存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th)這里Kf(Xn+th,Y(Xn+th)稱為平均斜率，龍格庫塔方法就是求得K的一種算法。利用這樣的原理，經(jīng)過復(fù)雜的數(shù)學(xué)推導(dǎo)（過于繁瑣省略），可以得出截斷誤差為O(h5)的四階龍格庫塔公式：K1f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K

3、4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)二、 Runge-Kutta的算法和流程圖在龍格-庫塔法中，四階龍格-庫塔法的局部截斷誤差約為0(h5),被廣泛應(yīng)用于解微分方程的初值問題。其算法公式為: 流程圖：(1)、四階龍格-庫塔方法流程圖: （2）、實(shí)例求解流程圖：三、 Runge-Kutta的Matlab實(shí)現(xiàn)function x,y=runge_kutta1(ufunc,y0,h,a,b)%參數(shù)表順序依次是微分方程組的函數(shù)名稱，初始值向量，步長，時間起點(diǎn)，時間終點(diǎn)（參數(shù)形式參考了ode45函數(shù)）n=floor(b-a)/h);%求步數(shù)x

4、(1)=a;%時間起點(diǎn)y(:,1)=y0;%賦初值，可以是向量，但是要注意維數(shù)for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii);k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龍格庫塔方法進(jìn)行數(shù)值求解end四、 Runge-Kutta的算例實(shí)現(xiàn)例:求解常微分方程 d(y)/d(x)=-2*y+2*x*x+2*x

5、 ，0x0.5 ，y(0)=1.編輯m文件Untitled3：fun=inline('-2*y+2*x*x+2*x');x,y=ode45(fun,0,0.5,1)>> Untitled3x = 0 0.0125 0.0250 0.0375 0.0500 0.0625 0.0750 0.0875 0.1000 0.1125 0.1250 0.1375 0.1500 0.1625 0.1750 0.1875 0.2000 0.2125 0.2250 0.2375 0.2500 0.2625 0.2750 0.2875 0.3000 0.3125 0.3250 0.3375 0.3500 0.3625 0.3750 0.3875 0.4000 0.4125 0.4250 0.4375 0.4500 0.4625 0.4750 0.4875 0.5000y = 1.0000 0.9755 0.9519 0.9291 0.9073 0.8864 0.8663 0.8471 0.8287 0.8112 0.7944 0.7785 0.7633 0.7489 0.7353 0.7224 0.7103 0.6989 0.6883 0.6783 0.6690 0.6605 0.6526 0.6454 0.6388 0.6329 0.6277 0.6