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

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

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

3、 求解一階常微分方程程序:yf = ft y) 譏 to) yowiki4LoveM 一階常微分方程定義函數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內|1C E.r1r-U>-JIjB2253J.B亠% ans =一階求解結果圖% (3*Ei(1) - 2*exp(1)/exp(1/t) - (3*Ei(1/t)/exp(1/t)2求解高階常微分方程關鍵是將高階轉為一階，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相關函數編輯本段回目錄ode23 , ode45 , ode113 , ode15s , ode23s , ode23t , ode23tbMatlab中龍格-庫塔(Runge-Kutta)方法原理及實現(xiàn)( 自己寫的，非直接調用 )龍格-庫塔(Ru nge-Kutta)方法是一種在工程上應用廣泛的高精度單步算法。由于 此算法精度高， 采取措施對誤差進行抑制， 所以其實現(xiàn)原理也較復雜。 該算法是 構建在數學支持的基礎之上的。 龍格庫塔方法的理論基礎來源于泰勒公式和使用 斜率近似表達微分，它在積分區(qū)間多預計算出幾個點的斜率， 然后進行加權平均， 用做下一點的依據，

7、從而構造出了精度更高的數值積分計算方法。 如果預先求兩 個點的斜率就是二階龍格庫塔法， 如果預先取四個點就是四階龍格庫塔法。 一階 常微分方程可以寫作： y'=f(x,y), 使用差分概念。(Yn+1-Yn)/h= f(Xn,Yn) 推出(近似等于，極限為 Yn') Yn+1=Yn+h*f(Xn,Yn)另外根據微分中值定理，存在 0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th)這里K = f(Xn+th,Y(Xn+th)稱為平均斜率，龍格庫塔方法就是求得K的一種算法。 利用這樣的原理，經過復雜的數學推導(過于繁瑣省略)，可以得出截斷誤差為 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) ; 所以，為了更好更準確地把握時間關系，應自己在理解龍格庫塔原理的基礎上， 編寫定步長的龍格庫塔函數， 經過學習其原理， 已經完成了一維的龍格庫塔函數。 仔細思考之后， 發(fā)現(xiàn)其實如果是需要解多個微分方程組， 可以想象成多個微分方 程并行進行求解， 時間，步長都是共同的， 首先把預定的初始值給每個微分方程 的第一步，然后每走一步，對多個微分

9、方程共同求解。想通之后發(fā)現(xiàn)，整個過程 其實很直觀，只是不停的逼近計算罷了。編寫的定步長的龍格庫塔計算函數：fun ctio n x,y=ru nge_kutta1(ufu nc,y0,h,a,b)%參數表順序依次是微分方程組的函數 名稱，初始值向量，步長，時間起點，時間終點(參數形式參考了ode45函數)n=floor(b-a)/h);% 求步數 x(1)=a;%時間起點y(:,1)=y0;%賦初值，可以是向量，但是要注意維數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;%按照龍格庫塔方法進行數值求解end 調用的子函數以及其調用語句： 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 和自編的龍格庫塔函數進行比較，調用如下： T,F = ode45(test_fun,0 15,1 1 3);subplot(121)plot(T,F)%Matlab自帶的ode45函數效果title(