ode45求解微分方程.doc

用MATLAB實(shí)現(xiàn)的四階龍格庫(kù)塔，并用來(lái)求解微分方程組function x,y=Runge_kutta(ufunc,y0,h,a,b)%參數(shù)表順序依次是微分方程組的函數(shù)名稱，初始值向量，步長(zhǎng)，時(shí)間起點(diǎn)，時(shí)間終點(diǎn)（參數(shù)形式參考了ode45函數(shù)）n=floor(b-a)/h); %求步數(shù),迭代次數(shù)%x=zeros(n+1,1); %n+1行1列的全零矩陣x(1)=a; %時(shí)間起點(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;%按照龍格庫(kù)塔方法進(jìn)行數(shù)值求解% y(:,ii+1)=roundn(y(:,ii+1),-8);endendclc;clear all;close all;%-初值和參數(shù)-a=16.923;b=20.653;r=0.51;G=1.41;k1=1/1.355;%k1=1/2.5;k2=1.25;k3=0;x0=0.001;0.1;0;0;%-定義方程-fun = (t,y) a*(y(2)-y(1)+G*y(1)-y(1)*( k1*(1-k2*(y(4)-k3)(-1/2) ); y(1)-y(2)+y(3); -b*y(2)-r*y(3); y(1); ;%-龍格庫(kù)塔解-t,x=Runge_Kutta(fun,x0,0.01,0,200);%測(cè)試時(shí)改變test_fun的函數(shù)維數(shù)，別忘記改變初始值的維數(shù)%-龍格庫(kù)塔解-v=x(1,3000:end);i=( (k1.*(1-k2.*(x(4,3000:end)-k3).(-1/2) ).*x(1,3000:end);x(1,:)=round(x(1,:)+13)*10+1);x(2,:)=round(x(2,:)+13)*10+1);x(3,:)=round(x(3,:)+13)*10+1);x(4,:)=round(x(4,:)+13)*10+1);v=fix(v+15)*10);i=fix(i+15)*10);figure(1)plot(t,x(1,:)%自編的龍格庫(kù)塔函數(shù)效果xlabel(x) ylabel(y);使用方法 T,Y = ode45(odefun,tspan,y0) odefun 是函數(shù)句柄,可以是函數(shù)文件名,匿名函數(shù)句柄或內(nèi)聯(lián)函數(shù)名 tspan 是區(qū)間 t0 tf 或者一系列散點(diǎn)t0,t1,.,tf y0 是初始值向量 T 返回列向量的時(shí)間點(diǎn) Y 返回對(duì)應(yīng)T的求解列向量 T,Y = ode45(odefun,tspan,y0,options) options 是求解參數(shù)設(shè)置,可以用odeset在計(jì)算前設(shè)定誤差,輸出參數(shù),事件等 T,Y,TE,YE,IE =ode45(odefun,tspan,y0,options) 每組(t,Y)之產(chǎn)生稱為事件函數(shù)。每次均會(huì)檢查是否函數(shù)等于零。并決定是否在零時(shí)終止運(yùn)算。這可以在函數(shù)中之特性上設(shè)定。例如以events 或events產(chǎn)生一函數(shù)。 value, isterminal,direction=events(t,y)其中，value(i)為函數(shù)之值，isterminal(i)=1時(shí)運(yùn)算在等于零時(shí)停止，=0時(shí)繼續(xù)；direction(i)=0時(shí)所有零時(shí)均需計(jì)算（默認(rèn)值）， +1在事件函數(shù)增加時(shí)等于零， -1在事件函數(shù)減少時(shí)等于零等狀況。此外，TE, YE, IE則分別為事件發(fā)生之時(shí)間，事件發(fā)生時(shí)之答案及事件函數(shù)消失時(shí)之指針i。 sol =ode45(odefun,t0 tf,y0.) sol 結(jié)構(gòu)體輸出結(jié)果 應(yīng)用舉例 1 求解一階常微分方程 程序: odefun=(t,y) (y+3*t)/t2; %定義函數(shù) tspan=1 4; %求解區(qū)間 y0=-2; %初值 t,y=ode45(odefun,tspan,y0); plot(t,y) %作圖 title(t2y=y+3t,y(1)=-2,1t4) legend(t2y=y+3t) xlabel(t) ylabel(y) % 精確解 % dsolve(t2*Dy=y+3*t,y(1)=-2) % ans = % (3*Ei(1) - 2*exp(1)/exp(1/t) - (3*Ei(1/t)/exp(1/t) 2 求解高階常微分方程 關(guān)鍵是將高階轉(zhuǎn)為一階,odefun的書(shū)寫(xiě). F(y,y,y.y(n-1),t)=0用變量替換,y1=y,y2=y.注意odefun方程定義為列向量 dxdy=y(1),y(2). 程序: function Testode45 tspan=3.9 4.0; %求解區(qū)間 y0=2 8; %初值 t,x=ode45(odefun,tspan,y0); plot(t,x(:,1),-o,t,x(:,2),-*) legend(y1,y2) title(y =-t*y + et*y +3sin2t) xlab