Matlab中解常微分方程的ode45

matlabode45解微分方程在用odesolver(ode45,ode15s,用法是：[t,y]=odesolver(odefun,tspan,y0);这里的odefunodefun用odesolver的时候传递参数呢？定义为global，这样做有一个弱点：针对系统不同，参数的表达与数用的时候进行传递，方法如下：的第四个参数一定是（odeset来定义option加进去，或者用使用空白矩阵（placeholder）。这样调用的时候格式就是：[t,y]=odesolver(odefun,tspan,y0,options,parameter1,parameter2);或者[t,y]=odesolver(odefun,tspan,y0,[],parameter1,parameter2);然后定义微分方程的时候也要有参数的地方：functiondydt=odefun(t,y,parameter1,parameter2)dydt=[eqn-1;eqn-2;…];就OK了。另一种用法：[T,Y]=ode45(@(t,y)rigid(t,y,Lambda,Lami,Cs,ionization,Ed),[0100],[PHI_0,E_0,M_0,Gma_0],options);Lambda,Lami,Cs,ionization,Edty是变量,[PHI_0,E_0,M_0,Gma_0],是初始值functiondy=rigid(t,y,Lambda,Lami,Cs,ionization,Ed)dy=zeros(4,1); %acolumnvectordy(1)=y(2);dy(2)=exp(y(1))-y(3)/y(4);dy(3)=Lambda*ionization/Cs*exp(y(1));dy(4)=-y(4)*Lambda/Lami-Lambda*ionization/Cs-y(2)/y(4)-Ed/y(4);~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~[t,H]=ode45('solitiontry1',tspan,h0,[],m1,epsinon)———————————————————————————————————————————————————————————functionsolfi=solitiontry1(t,H,flag,m1,epsinon)%m1=1.15;y=(1-epsinon)*(1-exp(H))+epsinon*m1^2*(1-sqrt(1+2*H*m1^-2))+m1^2*(1-sqrt(1-2*H*m1^-2));solfi=-sqrt(-2*y);flag参数。应该是针对options为空的那个参数的，但是在这里，如果options不为空的话，还不知道有什么区别。[原创[原创]Matlabode452010-09-2415:18ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kuttaode45表示采用四阶，五阶runge-kutta,(Δx)^3。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该,换用ode23来解.其他几个也是类似的用法使用方法[T,Y]=ode45(odefun,tspan,y0)odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan是区间[t0tf]或者一系列散点[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（默认值），+1-1YE,IE事件发生之时间事件发生时之答案事件函数消失时之指针i。sol=ode45(odefun,[t0sol=ode45(odefun,[t0tf],y0...)sol结构体输出结果应用举例1求解一阶常微分方程程序:odefun=@(t,y)(y+3*t)/t^2;%定义函数tspan=[14];%求解区间y0=-2;%初值[t,y]=ode45(odefun,tspan,y0);plot(t,y)%作图title('t^2y''=y+3t,y(1)=-2,1<t<4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')%精确解%dsolve('t^2*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). ]程序:functionTestode45tspan=[3.94.0];%求解区间y0=[28];%初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y''''=