Matlab中龍格-庫塔(Runge-Kutta)方法原理及實(shí)現(xiàn)

1、函數(shù)功能 編輯本段回目錄ode 是專門用于解微分方程的功能函數(shù)，他有 ode23,ode45,ode23s 等等，采用的是Runge-Kutta 算法。ode45表示采用四階，五階runge-kutta 單步算法,截斷誤差為(HX。 解決的是 Nonstiff( 非剛性 )的常微分方程 .是解決數(shù)值解問題的首選方法 ,若長時間沒結(jié)果 , 應(yīng)該就是剛性的 , 換用 ode23 來解 .使用方法 編輯本段回目錄T,Y = ode45(odefun,tspan,y0)odefun 是函數(shù)句柄 ,可以是函數(shù)文件名 ,匿名函數(shù)句柄或內(nèi)聯(lián)函數(shù)名tspan 是區(qū)間 t0 tf 或者一系列散點(diǎn) t0,t1,.

2、,tfy0是初始值向量T返回列向量的時間點(diǎn) www.iLoveMY返回對應(yīng) T 的求解列向量T,Y = ode45(odefun,tspan,y0,options)Simulink 與信號處理options 是求解參數(shù)設(shè)置 , 可以用 odeset 在計算前設(shè)定誤差 ,輸出參數(shù) ,事件等T,Y,TE,YE,IE =ode45(odefun,tspan,y0,options)在設(shè)置了事件參數(shù)后的對應(yīng)輸出 www.iLoveMTE事件發(fā)生時間 book.iLoveMYE事件解決時間IE事件消失時間sol =ode45(odefun,t0 tf,y0.)sol結(jié)構(gòu)體輸出結(jié)果應(yīng)用舉例 編輯本段回目錄1

3、 求解一階常微分方程程序:yf = ft y) 譏 to) yowiki4LoveM 一階常微分方程定義函數(shù)odefu n=(t,y) (y+3*t)/tA2;%tspa n=1 4;%求解區(qū)間y0=-2;%初值t,y=ode45(odefu n,tspa n,y 0);plot(t,y)%作圖title('tA2y”=y+3t,y(1)=-2,1<t<4')lege nd('tA2y"=y+3t')xlabel('t')Simuli nk與信號處理ylabel('y')%精確解 % dsolve('

4、tA2*Dy=y+3*t','y(1)=-2')« £Pypy+Jt 刖 TH內(nèi)|1C E.r1r-U>-JIjB2253J.B亠% ans =一階求解結(jié)果圖% (3*Ei(1) - 2*exp(1)/exp(1/t) - (3*Ei(1/t)/exp(1/t)2求解高階常微分方程關(guān)鍵是將高階轉(zhuǎn)為一階，odefun的書寫.F(y,y',y".y(n-1),t)=o用變量替換，y仁y,y2=y'.注意odefun方程定義為列向量dxdy=y(1),y(2)程序：function Testode45tspa n=3.9 4

5、.0; %求解區(qū)間y0=2 8;% 初值t,x=ode45(odefu n,tspa n,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')lege nd('y1','y2')title('y” ”=-t*y + eAt*y" +3si n2t') xlabel('t')ylabel('y')function y=odefu n( t,x) y=zeros(2,1); % 列向量 y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+

6、3*si n(2*t); end相關(guān)函數(shù)編輯本段回目錄ode23 , ode45 , ode113 , ode15s , ode23s , ode23t , ode23tbMatlab中龍格-庫塔(Runge-Kutta)方法原理及實(shí)現(xiàn)( 自己寫的，非直接調(diào)用 )龍格-庫塔(Ru nge-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ù)，

7、從而構(gòu)造出了精度更高的數(shù)值積分計算方法。 如果預(yù)先求兩 個點(diǎn)的斜率就是二階龍格庫塔法， 如果預(yù)先取四個點(diǎn)就是四階龍格庫塔法。 一階 常微分方程可以寫作： 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)這里K = f(Xn+th,Y(Xn+th)稱為平均斜率，龍格庫塔方法就是求得K的一種算法。 利用這樣的原理，經(jīng)過復(fù)雜的數(shù)學(xué)推導(dǎo)(過于繁瑣省略)，可以得出截斷誤差為 O(h

8、A5)的四階龍格庫塔公式：K1 = f(Xn，丫 n);K2=f(Xn+h/2,Yn+(h/2)*K1); K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6) ; 所以，為了更好更準(zhǔn)確地把握時間關(guān)系，應(yīng)自己在理解龍格庫塔原理的基礎(chǔ)上， 編寫定步長的龍格庫塔函數(shù)， 經(jīng)過學(xué)習(xí)其原理， 已經(jīng)完成了一維的龍格庫塔函數(shù)。 仔細(xì)思考之后， 發(fā)現(xiàn)其實(shí)如果是需要解多個微分方程組， 可以想象成多個微分方 程并行進(jìn)行求解， 時間，步長都是共同的， 首先把預(yù)定的初始值給每個微分方程 的第一步，然后每走一步，對多個微分

9、方程共同求解。想通之后發(fā)現(xiàn)，整個過程 其實(shí)很直觀，只是不停的逼近計算罷了。編寫的定步長的龍格庫塔計算函數(shù)：fun ctio n x,y=ru nge_kutta1(ufu nc,y0,h,a,b)%參數(shù)表順序依次是微分方程組的函數(shù) 名稱，初始值向量，步長，時間起點(diǎn)，時間終點(diǎn)(參數(shù)形式參考了ode45函數(shù))n=floor(b-a)/h);% 求步數(shù) x(1)=a;%時間起點(diǎn)y(:,1)=y0;%賦初值，可以是向量，但是要注意維數(shù)for ii=1:n x(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii);k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);

10、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 調(diào)用的子函數(shù)以及其調(diào)用語句： function dy=test_fun(x,y)dy = zeros(3,1);%初始化列向量dy(1) = y(2) * y(3);dy(2) = -y(1) + y(3);dy(3) = -0.51 * y(1) * y(2);對該微分方程組用 ode45 和自編的龍格庫塔函數(shù)進(jìn)行比較，調(diào)用如下： T,F = ode45(test_fun,0 15,1 1 3);subplot(121)plot(T,F)%Matlab自帶的ode45函數(shù)效果title(