matlab變參量微分方程參數(shù)識(shí)別(共13頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上1 變參數(shù)微分方程數(shù)值求解例子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 匿名函數(shù)法f=(t,y,u,v) u+2

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)-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 修改加上時(shí)間tt(顯示所有計(jì)算值)clearu=1;5;15;20;25;v=6;12;18;24;30;tspan=0:1:4;y0=0 2;tt =;yy=;for i=1:l

3、ength(tspan)-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);%所有的計(jì)算數(shù)值。2.3 同過(guò)差值可以調(diào)節(jié)精度。如果u,v隨著t是時(shí)刻變化的，但是通過(guò)測(cè)試手段只能測(cè)得某一時(shí)刻的u,v. clearglobal yyt1=0:0.1:4;%如果u,v隨著t是時(shí)刻變化的，可以通過(guò)此數(shù)值來(lái)調(diào)節(jié)精度tspan=0:1:4;u=1;5;15;20;25; u1=spline(tspan,u,t1);v=6;12;18;24;30;v1= spline(tspa

4、n,v,t1);y0=0 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對(duì)一個(gè)常微分方程進(jìn)行參數(shù)回歸問(wèn)題如下：已知實(shí)驗(yàn)數(shù)據(jù)x,y，并且x,y的關(guān)系滿(mǎn)足以下常微分方程Dy/dx=-k*(y-y0)*y2其中 k是需要回歸的參數(shù)，y0是一個(gè)常數(shù)，通常等于y向量中的最后一個(gè)數(shù)值。要求：1.通過(guò)lsqcurvefit或者lsqnonlin回歸出系數(shù)k2.畫(huà)出模型預(yù)測(cè)值和實(shí)驗(yàn)值的對(duì)比圖，模型預(yù)測(cè)值可以通過(guò)得到常微分方程

5、數(shù)值解后三次樣條spline插值得到。我已經(jīng)寫(xiě)好的程序如下，里面有錯(cuò)誤，我自己找不出來(lái)，請(qǐng)高手幫幫忙，謝謝啊可以加我的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,&

6、#39;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);

7、y=yc;function dydt = modeleqs(t,y,beta)global Jedydt =-beta*(y-Je)*y.2;=數(shù)據(jù)如下：X    Y0        1176.1        637.2        326.3        265.4        220.5        200.7 

8、      181.9        170.11        162.14        151.17        150.20        146.25        139.30        137.35        135.40&#

9、160;       131.45        131.53        126.60        125.70        123.80        121.90        120.100        116.110       

10、117.120        114.130        113.140        113.150        112.3變參數(shù)微分方程系數(shù)回歸您好斑竹，我的最終目的是要識(shí)別k參數(shù)，我的微分方程模型為y'=k(y-y0)*y2+v ;(我的參數(shù)識(shí)別方法是參考一個(gè)例子，后面給貼出來(lái)。但是這個(gè)例子是不帶v輸入變量的）1 我前面一段程序先假設(shè)k=3,然后利用變參量數(shù)值解法求解了（仿真）y值；(因?yàn)檫@個(gè)式子有個(gè)帶外部輸入的參數(shù)v

11、，沒(méi)法求得解析解，所以只好用數(shù)值微分方程方法求出ydata；2 我后面一段程序就是利用前邊得到的t值和y值去識(shí)別參數(shù)k；3.1 生成數(shù)據(jù)xxdata為2 的x數(shù)據(jù)。function dydt=modelsq(t,y,k,v)dydt =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;%計(jì)算出500個(gè)數(shù)據(jù)；t1=min(x):(max(x)-min

12、(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),y0,k,v(i); y0=y(end,: );yy=yy;y0;endxydata=t1',yy,v'save('xydata.txt','xydata','-ascii');%把tt和yy存為txt文件。3.2 參數(shù)識(shí)別方法1 lsqcurvefitfunction dydt=modelsq(t,y,k,v

13、)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) 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

14、('TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',300, 'Algorithm','trust-region-reflective', '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;

15、load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);F =(k)num_value(k,xdata)-ydata;ms=MultiStart('Display','iter', 'UseParallel','always');matlabpool open 2;opts=optimset('Algorithm','trust-region-reflective');problem=createOptimProblem(&

16、#39;lsqnonlin','x0',-4,2,'objective',F,'lb', -6 1,'ub',-1 5,'options',opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)matlabpool closeopts=optimset('Algorithm','trust-region-reflective','MaxTime',200，'Maxiter',40

17、0);方法3 lsqcurvefit 和multistartglobal v;load xydata.txt;xdata=xydata(:,1);ydata=xydata(:,2);v=xydata(:,3);ms=MultiStart('Display','iter','UseParallel','always');opts=optimset('Algorithm','trust-region-reflective');matlabpool openproblem=createOptimProbl

18、em('lsqcurvefit','x0',-4,2, 'objective',num_value,'xdata',xdata,'ydata',ydata,'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=xy

19、data(:,2);v=xydata(:,3);F=(k)norm(num_value(k,xdata)-ydata);gs=GlobalSearch('Display','iter');matlabpool openopts=optimset('Algorithm','interior-point', 'UseParallel','always');problem=createOptimProblem('fmincon','x0',-4,2,'objecti

20、ve',F,'lb',-6 1,'ub',-1 5,'options',opts); xming,fming,flagg,outptg,manyminsg = run(gs,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= saop

21、timset('Display','iter');x,fval,exitFlag,output = simulannealbnd(F,X0,lb,ub,options)time=toc;二、關(guān)于求解變參數(shù)微分方程看到有不少人問(wèn)過(guò)如何處理"變參數(shù)微分方程", 所以抽了一點(diǎn)時(shí)間,寫(xiě)了一個(gè)例子,希望能幫助到那些需要求解此類(lèi)問(wèn)題的人。%=%clear allfun=inline('y(2);sin(w*t)-2*y(1)-3*y(2)','t','y','flag','w'

22、;); tsp=0 10;y0=1 1;xlim(tsp) for w=1:10   t,y=ode45(fun,tsp,y0,w);    plot(t,y)   title('w = ',num2str(w);    pause(0.5); end三、關(guān)于求解變上限積分問(wèn)題經(jīng)?？吹接腥嗽?xún)問(wèn)如何求解“變上限積分問(wèn)題”，現(xiàn)舉一個(gè)例子，說(shuō)明該如何求解之：希望以后碰到類(lèi)似問(wèn)題時(shí)，能夠自己動(dòng)手解決%=%clear allwarning offx=linspace(0,150);for k=1:length

23、(x)    xx=x(k);    fun=inline('exp(-lam.2)');    w(k)=0.62*2/sqrt(pi)*quadl(fun,0,sqrt(pi)*xx/150); endplot(x,w)4 變系數(shù)數(shù)值積分問(wèn)題。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);%

24、這個(gè)實(shí)際上是四缸機(jī)的扭矩要求Ti在0,4*pi,上的數(shù)值積分，而且每隔一固定角度間隔，要從外部輸入p1,p2,p3,p4;function y=calTi(o,p1,p2,p3,p4) R =52.5e-3;%曲柄半徑L=184e-3;%連桿長(zhǎng)度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

25、 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(tspan)-1z=quadl(calTi,tspan(i),tspan(i+1),p1(i),p2(i), p3(i), p4(i);Ti1(i, :)=z end5.1關(guān)于帶參數(shù)的積分問(wèn)題例如以下問(wèn)題: 函數(shù)為 y=sin(k.*x).*x.2，對(duì)x積分，積分區(qū)域?yàn)椤?，5】，目的是要畫(huà) k 和 y 的圖形.方法一（用循環(huán)）：1. f=(k) q

26、uad(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. plot(kk,y)復(fù) 制代碼方法二（不用循環(huán)）：1. plot(linspace(0,5),arrayfun(k) quad(x)  sin(k.*x).*x.2,0,5),linspace(0,5)復(fù)制代碼振動(dòng)方程求解到有不少人問(wèn)過(guò)二階動(dòng)力微分方程的求解問(wèn)題,現(xiàn)舉一個(gè)簡(jiǎn)單的例子, 其余的情形希望讀者能舉一反三, 自己多思考.%=%clear al

27、ln=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=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)',.    

28、;      't','x','flag','n','M','K','F','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函數(shù)的使用自己的實(shí)驗(yàn)；1

29、t=0:0.1:5; T=arrayfun(x) x.2,t); 0-5的挨個(gè)平方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 求f(k)=int(sin(k*x)*x2,0,5,x),k取0,5上的不同值；方法1 匿名函數(shù)tic;f=(k)quadl(x)sin(k.*x).*x.2,0,5);%相當(dāng)于多重匿名函數(shù)呵呵。kk=linspace(0,50;%默認(rèn)分為100等分y2=zeros(size(kk);for ii=1:length(kk)y2(ii)=f(kk(ii);endtime=

30、toc;方法2：向量化方法+匿名函數(shù)tic ;y=arrayfun(k) quadl(x) sin(k*x).*x.2,0,5),linspace(0,5);%相當(dāng)于求一百組循環(huán)代入k。time=toc;例子1 匿名函數(shù)數(shù)值微分方程f=(t,x)x(2); -t*x(1)+exp(t)*x(2)+3*sin(2*t);tspan=3.9 4.0; %求解區(qū)間 y0=2 8; %初值t,x=ode45(f,tspan,y0); plot(t,x(:,1),'-o',t,x(:,2),'-*') ;legend('y1','y2');

31、title('y'' ''=-t*y + et*y'' +3sin2t') ;xlabel('t') ;ylabel('y');例子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 多變量匿名函數(shù)a=1;b=2;g=(x,y)a*x+y.b;g1:5,1:5);3.1 f=(x,y)x.2+y.2; %注意需要點(diǎn)(.)運(yùn)算a=1:1:10;b=10:-1:1;f(a,

32、b)3.2匿名函數(shù)的表達(dá)式中也可以有參數(shù)的傳遞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 多重匿名函數(shù)f=(x,y)(a) x2+y+a;f1=f(2,3)f1 =(a)x2+y+af2=f1(4)f2 =854.1 匿名函數(shù)和m函數(shù)的聯(lián)合function e=myfun(x,y)e=2.*exp(4.*x)+8.*sin(x)-y;feval(x,y)myfun(x,y),2,3)4.2使用匿名函數(shù)實(shí)現(xiàn)符號(hào)函數(shù)的賦值運(yùn)算>> syms x; %定義符號(hào)變量>> z=2*x3+4*x+5; %定義表達(dá)式&g

33、t;> z1=diff(z,2) %求z的2階導(dǎo)數(shù)的表達(dá)式z1 =12*xz2=eval('(x)' vectorize(z1); %vectorize函數(shù)的功能是使內(nèi)聯(lián)函數(shù)適合數(shù)組運(yùn)算的法則>> z2(3)ans=36例子3 匿名函數(shù)的求極值f=(x)exp(x)+x2+x(sqrt(x)-100;format long ;x0=fzero(f,3);x0;f(x0);3.1求aa在不同值時(shí)下列方程的根f=(a)(x)exp(x)+xa+x(sqrt(x)-100;aa=0:0.01;2;arrayfun(a)fzero(f(a),4),aa);%4為搜索的

34、初值。例子4 匿名函數(shù)在現(xiàn)實(shí)表示隱函數(shù)y=(x)fzero(y)(exp(y)+xy)(1/y)-x2*y,1);這樣對(duì)于任意的x，只要調(diào)用y(x),就能得到相應(yīng)的y值,如y(1)這時(shí)y只接受標(biāo)量，利用arrayfun函數(shù)Y=(x)arrayfun(xxx)fzero(y)(exp(y)+xy)(1/y)-x2*y,1),x);Y(1:10);例子4 利用匿名函數(shù)的數(shù)值積分quad2d(x,y)x.*y,1,2,(x)sin(x),(x)cos(x),AbsTol,1e-12);4.1可以接受向量的表達(dá)式quadl(x)arrayfun(xx)quadl(y)xx*y,sin(xx),cos(

35、xx),x),1,2);5 matlab 改進(jìn)的隨機(jī)行走法%x為初始值，lamda為步長(zhǎng)，epsilon 為控制迭代參數(shù)，即lamda減小到epsilon時(shí)迭代終止；mx為n實(shí)驗(yàn)得到的最優(yōu)解，minf為相應(yīng)的最優(yōu)質(zhì)。function 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

36、(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);結(jié)果mx=50;minf=-1.2492;例6 如求以下這個(gè)函數(shù)f（x(1)-2).2+(x(2)-3