ode45求解微分方程.doc

用MATLAB实现的四阶龙格库塔，并用来求解微分方程组function x,y=Runge_kutta(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor(b-a)/h); %求步数,迭代次数%x=zeros(n+1,1); %n+1行1列的全零矩阵x(1)=a; %时间起点y(:,1)=y0; %赋初值，可以是向量，但是要注意维数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;%按照龙格库塔方法进行数值求解% y(:,ii+1)=roundn(y(:,ii+1),-8);endendclc;clear all;close all;%-初值和参数-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); ;%-龙格库塔解-t,x=Runge_Kutta(fun,x0,0.01,0,200);%测试时改变test_fun的函数维数，别忘记改变初始值的维数%-龙格库塔解-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,:)%自编的龙格库塔函数效果xlabel(x) ylabel(y);使用方法 T,Y = ode45(odefun,tspan,y0) odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 tspan 是区间 t0 tf 或者一系列散点t0,t1,.,tf y0 是初始值向量 T 返回列向量的时间点 Y 返回对应T的求解列向量 T,Y = ode45(odefun,tspan,y0,options) options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 T,Y,TE,YE,IE =ode45(odefun,tspan,y0,options) 每组(t,Y)之产生称为事件函数。每次均会检查是否函数等于零。并决定是否在零时终止运算。这可以在函数中之特性上设定。例如以events 或events产生一函数。 value, isterminal,direction=events(t,y)其中，value(i)为函数之值，isterminal(i)=1时运算在等于零时停止，=0时继续；direction(i)=0时所有零时均需计算（默认值）， +1在事件函数增加时等于零， -1在事件函数减少时等于零等状况。此外，TE, YE, IE则分别为事件发生之时间，事件发生时之答案及事件函数消失时之指针i。 sol =ode45(odefun,t0 tf,y0.) sol 结构体输出结果 应用举例 1 求解一阶常微分方程 程序: odefun=(t,y) (y+3*t)/t2; %定义函数 tspan=1 4; %求解区间 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 求解高阶常微分方程 关键是将高阶转为一阶,odefun的书写. 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; %求解区间 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