matlab变参量微分方程处理.doc

1 变参数微分方程数值求解 例子function dydt=fun(t,y,u,v) r=u+2;s=v-2; dydt=[r+y(2); s*y(1)-2*s*y(2)]; u=[1;5;15;20;25];v=[6;12;18;24;30]; tspan=0:1:4;y0=[0 2]; yy=y0;for i=1:length(tspan)-1[t,y]=ode45(@fun,[tspan(i),tspan(i+1)],y0,[],u(i),v(i));y0=y(end,: );yy=[yy;y0];endplot(tspan,yy,-o) 2.1 匿名函数法 f=@(t,y,u,v) [u+2+y(2); (v-2)*y(1)-2*(v-2)*y(2)]u=[1;5;15;20;25];v=[6;12;18;24;30]; tspan=0:1:4;y0=[0 2]; yy=y0;for i=1:length(tspan)-1[t,y]=ode45(f,[tspan(i),tspan(i+1)],y0,[],u(i),v(i));y0=y(end,: );yy=[yy;y0];endplot(tspan,yy,-o) 2.2 修改加上时间tt(显示所有计算值) clearu=[1;5;15;20;25];v=[6;12;18;24;30]; tspan=0:1:4;y0=[0 2]; tt =[];yy=[];for i=1:length(tspan)-1[t,y]=ode45(@fun,[tspan(i),tspan(i+1)],y0,[],u(i),v(i));y0=y(end,: );tt=[tt;t]; yy=[yy;y];end plot(tt,yy);%所有的计算数值。 2.3 同过差值可以调节精度。如果u,v随着t是时刻变化的，但是通过测试手段只能测得某一时刻的u,v. clear global yy t1=0:0.1:4;%如果u,v随着t是时刻变化的，可以通过此数值来调节精度 tspan=0:1:4;u=[1;5;15;20;25]; u1=spline(tspan,u,t1);v=[6;12;18;24;30];v1= spline(tspan,v,t1);y0=[0 2]; yy=y0;for i=1:length(t1)-1[t,y]=ode45(@fun,[t1(i),t1(i+1)],y0,[],u1(i),v1(i));y0=y(end,: ); yy=[yy;y0];end plot(t1,yy)2 适用matlab对一个常微分方程进行参数回归 问题如下：已知实验数据x,y，并且x,y的关系满足以下常微分方程Dy/dx=-k*(y-y0)*y^2其中 k是需要回归的参数，y0是一个常数，通常等于y向量中的最后一个数值。要求：1.通过lsqcurvefit或者lsqnonlin回归出系数k2.画出模型预测值和实验值的对比图，模型预测值可以通过得到常微分方程数值解后三次样条spline插值得到。我已经写好的程序如下，里面有错误，我自己找不出来，请高手帮帮忙，谢谢啊可以加我的QQ交流=======================================function odetestclc;clear;global Je J0data=xlsread(flux.xls);xdata=data(:,1);ydata=data(:,2);beta0=0.1;Je=ydata(end);J0=ydata(1);options=optimset(TolFun,1e-20,TolX,1e-20,MaxFunEvals,100,Algorithm,trust-region-reflective); beta=lsqcurvefit(@cakefun,beta0,xdata,ydata,[],[],options);Jc=cakefun(beta,xdata);plot(xdata,ydata,o,xdata,Jc);function y=cakefun(beta,x)global J0tspan=[0 max(x)];[m,n]=size(x);[tt yy] = ode23s(@modeleqs,tspan,J0,[],beta); yc=spline(tt,yy,x);y=yc;function dydt = modeleqs(t,y,beta)global Jedydt =-beta*(y-Je)*y.^2;=====================