matlab变参量微分方程参数识别

1、1 变参数微分方程数值求解例子2 求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)-1t,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)*

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)-1t,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)-

3、1t,y=ode45(fun,tspan(i),tspan(i+1),y0,u(i),v(i);y0=y(end,: );tt=tt;t;yy=yy;y;endplot(tt,yy);%所有的计算数值。2.3 同过差值可以调节精度。如果u,v随着t是时刻变化的，但是通过测试手段只能测得某一时刻的u,v. clearglobal yyt1=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

4、2;yy=y0;for i=1:length(t1)-1t,y=ode45(fun,t1(i),t1(i+1),y0,u1(i),v1(i);y0=y(end,: );yy=yy;y0;endplot(t1,yy)2 适用matlab对一个常微分方程进行参数回归问题如下：已知实验数据x,y，并且x,y的关系满足以下常微分方程Dy/dx=-k*(y-y0)*y2其中 k是需要回归的参数，y0是一个常数，通常等于y向量中的最后一个数值。要求：1.通过lsqcurvefit或者lsqnonlin回归出系数k2.画出模型预测值和实验值的对比图，模型预测值可以通过得到常微分方程数值解后三次样条splin

5、e插值得到。我已经写好的程序如下，里面有错误，我自己找不出来，请高手帮帮忙，谢谢啊可以加我的QQ交流：40231185=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,'A

6、lgorithm','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;

7、function dydt = modeleqs(t,y,beta)global Jedydt =-beta*(y-Je)*y.2;=数据如下：X    Y0        1176.9181151        637.42257272        326.12181033        265.76315284        220.6660835   

8、;     200.79882657        181.62981219        170.163443611        162.468402414        151.632275917        150.381132820        146.824206925       

9、 139.873516430        137.586118635        135.109397740        131.719599445        131.63195153        126.415986560        125.021992670        123.04196780&

10、#160;       121.734474190        120.8522371100        116.8166294110        117.2211892120        114.8271487130        113.2291363140        113.2365511150  &

11、#160;     112.96558663变参数微分方程系数回归您好斑竹，我的最终目的是要识别k参数，我的微分方程模型为y'=k(y-y0)*y2+v ;(我的参数识别方法是参考一个例子，后面给贴出来。但是这个例子是不带v输入变量的）1 我前面一段程序先假设k=3,然后利用变参量数值解法求解了（仿真）y值；(因为这个式子有个带外部输入的参数v，没法求得解析解，所以只好用数值微分方程方法求出ydata；2 我后面一段程序就是利用前边得到的t值和y值去识别参数k；3.1 生成数据xxdata为2 的x数据。function dydt=modelsq(t,y,k,v)d

12、ydt =k(1)*(y-6).k(2)+v;clear;k=-3.548 3.236;y0=6;x=0 1 2 3 4 5 7 9 11 14 17 20 25 30 35 40 45 53 60 70 80 90 100 110 120 130 140 150'%给出x值fraction=500;%计算出500个数据；t1=min(x):(max(x)-min(x)/fraction:max(x);v1=1:length(x);v=spline(x,v1,t1);yy=y0;for i=1:length(t1)-1t y=ode23s(modelsq,t1(i) t1(i+1),y

13、0,k,v(i); y0=y(end,: );yy=yy;y0;endxydata=t1',yy,v'save('xydata.txt','xydata','-ascii');%把tt和yy存为txt文件。3.2 参数识别方法1 lsqcurvefitfunction dydt=modelsq(t,y,k,v)dydt =k(1)*(y-6).k(2)+v;function yi=num_value(k,x)global v;y0=6;yy=y0;for i=1:length(x)-1t y=ode23s(modelsq,x(i)

14、 x(i+1),y0, ,k,v(i); y0=y(end,: );yy=yy;y0;endyi=yy;clear;global v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);k0=4000,2000' %初值也可以改为k0=-6,2'options=optimset('TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',300, 'Algorithm','trust-region-refle

15、ctive', 'display', 'iter');lb=-13 1;ub=-1 5;kp=lsqcurvefit(num_value,k0,xdata,ydata,lb,ub,options);%kp =-3.5400,3.23000;方法2 multistart和lsqnonglin, 全局算法（速度比较慢）global v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);F =(k)num_value(k,xdata)-ydata;ms=MultiStart(&#

16、39;Display','iter', 'UseParallel','always');matlabpool open 2;opts=optimset('Algorithm','trust-region-reflective');problem=createOptimProblem('lsqnonlin','x0',-4,2,'objective',F,'lb', -6 1,'ub',-1 5,'options'

17、,opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)matlabpool closeopts=optimset('Algorithm','trust-region-reflective','MaxTime',200，'Maxiter',400);方法3 lsqcurvefit 和multistartglobal v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);ms=Multi

18、Start('Display','iter','UseParallel','always');opts=optimset('Algorithm','trust-region-reflective');matlabpool openproblem=createOptimProblem('lsqcurvefit','x0',-4,2, 'objective',num_value,'xdata',xdata,'ydata',y

19、data,'lb',-6,1,'ub',-1,5,'options',opts); xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)方法4 globalsearchglobal v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);F=(k)norm(num_value(k,xdata)-ydata);gs=GlobalSearch('Display','iter');mat

20、labpool openopts=optimset('Algorithm','interior-point', 'UseParallel','always');problem=createOptimProblem('fmincon','x0',-4,2,'objective',F,'lb',-6 1,'ub',-1 5,'options',opts); xming,fming,flagg,outptg,manyminsg = run(g

21、s,problem)方法5 模拟退火tic;global v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);F=(k)norm(num_value(k,xdata)-ydata);X0 =-4,2;% Starting pointlb=-6 1;ub=-1 5;options= saoptimset('Display','iter');x,fval,exitFlag,output = simulannealbnd(F,X0,lb,ub,options)time=toc;二、关于

22、求解变参数微分方程看到有不少人问过如何处理"变参数微分方程", 所以抽了一点时间,写了一个例子,希望能帮助到那些需要求解此类问题的人。%=%clear allfun=inline('y(2);sin(w*t)-2*y(1)-3*y(2)','t','y','flag','w'); tsp=0 10;y0=1 1;xlim(tsp) for w=1:10   t,y=ode45(fun,tsp,y0,w);    plot(t,y)   

23、title('w = ',num2str(w);    pause(0.5); end三、关于求解变上限积分问题经常看到有人询问如何求解“变上限积分问题”，现举一个例子，说明该如何求解之：希望以后碰到类似问题时，能够自己动手解决%=%clear allwarning offx=linspace(0,150);for k=1:length(x)    xx=x(k);    fun=inline('exp(-lam.2)');    w(k)=0.62*2/sqrt(pi)*quadl

24、(fun,0,sqrt(pi)*xx/150); endplot(x,w)4 变系数数值积分问题。Ti=Ap*R*(p1+p3).*(sin(o)+(R*sin(2*o)./(2*L*(1-(R2*sin(o).2)/L2).(1/2)-Ap*R*(p2+p4).*(sin(o)-(R*sin(2*o)./(2*L*(1-(R2*sin(o).2)/L2).(1/2);%这个实际上是四缸机的扭矩要求Ti在0,4*pi,上的数值积分，而且每隔一固定角度间隔，要从外部输入p1,p2,p3,p4;function y=calTi(o,p1,p2,p3,p4) R =52.5e-3;%曲柄半径L=18

25、4e-3;%连杆长度Ap= 0.0071;%活塞面积y=Ap*R*(p1+p3).*(sin(o)+(R*sin(2*o)./(2*L*(1-(R2*sin(o).2)/L2).(1/2)-Ap*R*(p2+p4).*(sin(o)-(R*sin(2*o)./(2*L*(1-(R2*sin(o).2)/L2).(1/2);%matlab空间代码clc;clear;load ptwo.txt;p1=ptwo(:,1); p2=ptwo(:,2); p3=ptwo(:,3); p4=ptwo(:,4);tspan=0:4*pi/length(p1):4*pi;Ti1=;for i=1:length

26、(tspan)-1z=quadl(calTi,tspan(i),tspan(i+1),p1(i),p2(i), p3(i), p4(i);Ti1(i, :)=z end5.1关于带参数的积分问题例如以下问题: 函数为 y=sin(k.*x).*x.2，对x积分，积分区域为【1，5】，目的是要画 k 和 y 的图形.方法一（用循环）：1. f=(k) quad(x)sin(k.*x).*x.2,0,5)2. kk=linspace(0,5);3. y=zeros(size(kk);4. for ii=1:length(kk)5. y(ii)=f(kk(ii);6. end7.

27、 plot(kk,y)复 制代码方法二（不用循环）：1. plot(linspace(0,5),arrayfun(k) quad(x)  sin(k.*x).*x.2,0,5),linspace(0,5)复制代码振动方程求解到有不少人问过二阶动力微分方程的求解问题,现举一个简单的例子, 其余的情形希望读者能举一反三, 自己多思考.%=%clear alln=3;F=25;24;20;m1=31.2; m2=31.2;  m3=31.2;k1=67.51;k2=67.51; k3=67.51;c1=0.01; c2=0.01;  c3=

28、0.01;M=m1,0,0;0,m2,0;0,0,m3;B=c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3;K=k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3;DL=inline('x(n+1:end,1); inv(M)*(F-B*x(1:n,1)-K*x(1:n,1)',.          't','x','flag','n','M','K','F','

29、;B');options = odeset('RelTol',1e-4,'AbsTol',1e-4 1e-4 1e-5);t,x=ode45(DL,0,3,rand(n,1),options,n,M,K,F,B);  plot(t,x(:,1:n)MatLab向量化编程:arrayfun函数的使用自己的实验；1 t=0:0.1:5; T=arrayfun(x) x.2,t); 0-5的挨个平方2 t=0:0.1:5; f=(x)x.2; f(t);同上3 t=0:0.1:5; f=(x)x.2; y=arrayfun(f,t);同上1

30、 求f(k)=int(sin(k*x)*x2,0,5,x),k取0,5上的不同值；方法1 匿名函数tic;f=(k)quadl(x)sin(k.*x).*x.2,0,5);%相当于多重匿名函数呵呵。kk=linspace(0,50;%默认分为100等分y2=zeros(size(kk);for ii=1:length(kk)y2(ii)=f(kk(ii);endtime=toc;方法2：向量化方法+匿名函数tic ;y=arrayfun(k) quadl(x) sin(k*x).*x.2,0,5),linspace(0,5);%相当于求一百组循环代入k。time=toc;例子1 匿名函数数值微

31、分方程f=(t,x)x(2); -t*x(1)+exp(t)*x(2)+3*sin(2*t);tspan=3.9 4.0; %求解区间 y0=2 8; %初值t,x=ode45(f,tspan,y0); plot(t,x(:,1),'-o',t,x(:,2),'-*') ;legend('y1','y2');title('y'' ''=-t*y + et*y'' +3sin2t') ;xlabel('t') ;ylabel('y');例

32、子2 求值f=(a,b,x) a*x+b; a=1:0.1:2;b=2:0.1:3;x=1:10;X,Y,Z=meshgrid(a,b,x);V=arrayfun(f,X,Y,Z);3 多变量匿名函数a=1;b=2;g=(x,y)a*x+y.b;g1:5,1:5);3.1 f=(x,y)x.2+y.2; %注意需要点(.)运算a=1:1:10;b=10:-1:1;f(a,b)3.2匿名函数的表达式中也可以有参数的传递a=1:5;b=5:-1:1;c=0.1:0.1:0.5;f=(x,y)x.2+y.2+c;f(a,b)4 多重匿名函数f=(x,y)(a) x2+y+a;f1=f(2,3)f1

33、=(a)x2+y+af2=f1(4)f2 =854.1 匿名函数和m函数的联合function e=myfun(x,y)e=2.*exp(4.*x)+8.*sin(x)-y;feval(x,y)myfun(x,y),2,3)4.2使用匿名函数实现符号函数的赋值运算>> syms x; %定义符号变量>> z=2*x3+4*x+5; %定义表达式>> z1=diff(z,2) %求z的2阶导数的表达式z1 =12*xz2=eval('(x)' vectorize(z1); %vectorize函数的功能是使内联函数适合数组运算的法则>&g

34、t; z2(3)ans=36例子3 匿名函数的求极值f=(x)exp(x)+x2+x(sqrt(x)-100;format long ;x0=fzero(f,3);x0;f(x0);3.1求aa在不同值时下列方程的根f=(a)(x)exp(x)+xa+x(sqrt(x)-100;aa=0:0.01;2;arrayfun(a)fzero(f(a),4),aa);%4为搜索的初值。例子4 匿名函数在现实表示隐函数y=(x)fzero(y)(exp(y)+xy)(1/y)-x2*y,1);这样对于任意的x，只要调用y(x),就能得到相应的y值,如y(1)这时y只接受标量，利用arrayfun函数Y=

35、(x)arrayfun(xxx)fzero(y)(exp(y)+xy)(1/y)-x2*y,1),x);Y(1:10);例子4 利用匿名函数的数值积分quad2d(x,y)x.*y,1,2,(x)sin(x),(x)cos(x),AbsTol,1e-12);4.1可以接受向量的表达式quadl(x)arrayfun(xx)quadl(y)xx*y,sin(xx),cos(xx),x),1,2);5 matlab 改进的随机行走法%x为初始值，lamda为步长，epsilon 为控制迭代参数，即lamda减小到epsilon时迭代终止；mx为n实验得到的最优解，minf为相应的最优质。funct

36、ion mx,minf=randwalk(f,x,lamda,epsilon,N,n)F=zeros(n,1);X=zeros(n,2);f1=f(x(1),x(2);while(lamda>=epsilon)k=1;while(k<=N)u=unifrnd(-1,1,n,2);for ii=1:nX(ii,:)=x+lamda*(u(ii,:)/norm(u(ii,:);F(ii)=f(X(ii,1),X(ii,2);endf11,kk=min(F);if f11<f1f1=f11;x=X(kk,:);k=1;else k=k+1end endlamda=lamda/2;endmx=X(kk,:);minf=f1;f=(x,y)-sin(sqrt(x-50).2+(y-50).2)+exp(1)./(sqrt(x-50).2+(y-50).2)+exp(1)-1;mx,minf=randwalk(f, 0,0, 10, 0.00001, 1000, 10);结果mx=50;minf=-1.2492;例6 如