matlab变参量微分方程处理.docx
- 文档编号:4404055
- 上传时间:2022-12-01
- 格式:DOCX
- 页数:16
- 大小:35.49KB
matlab变参量微分方程处理.docx
《matlab变参量微分方程处理.docx》由会员分享,可在线阅读,更多相关《matlab变参量微分方程处理.docx(16页珍藏版)》请在冰豆网上搜索。
matlab变参量微分方程处理
1变参数微分方程数值求解
例子2求
functiondydt=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=[02];
yy=y0;
fori=1:
length(tspan)-1
[t,y]=ode45(@fun,[tspan(i),tspan(i+1)],y0,[],u(i),v(i));
y0=y(end,:
);
yy=[yy;y0];
end
plot(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=[02];
yy=y0;
fori=1:
length(tspan)-1
[t,y]=ode45(f,[tspan(i),tspan(i+1)],y0,[],u(i),v(i));
y0=y(end,:
);
yy=[yy;y0];
end
plot(tspan,yy,'-o')
2.2修改加上时间tt(显示所有计算值)
clear
u=[1;5;15;20;25];
v=[6;12;18;24;30];
tspan=0:
1:
4;
y0=[02];
tt=[];yy=[];
fori=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
globalyy
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=[02];
yy=y0;
fori=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回归出系数k
2.画出模型预测值和实验值的对比图,模型预测值可以通过得到常微分方程数值解后三次样条spline插值得到。
我已经写好的程序如下,里面有错误,我自己找不出来,请高手帮帮忙,谢谢啊
可以加我的QQ交流:
40231185
=======================================
functionodetest
clc;clear;
globalJeJ0
data=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);
functiony=cakefun(beta,x)
globalJ0
tspan=[0max(x)];
[m,n]=size(x);
[ttyy]=ode23s(@modeleqs,tspan,J0,[],beta);
yc=spline(tt',yy',x);
y=yc;
functiondydt=modeleqs(t,y,beta)
globalJe
dydt=-beta*(y-Je)*y.^2;
===========================
数据如下:
X Y
0 1176.918115
1 637.4225727
2 326.1218103
3 265.7631528
4 220.666083
5 200.7988265
7 181.6298121
9 170.1634436
11 162.4684024
14 151.6322759
17 150.3811328
20 146.8242069
25 139.8735164
30 137.5861186
35 135.1093977
40 131.7195994
45 131.631951
53 126.4159865
60 125.0219926
70 123.041967
80 121.7344741
90 120.8522371
100 116.8166294
110 117.2211892
120 114.8271487
130 113.2291363
140 113.2365511
150 112.9655866
3变参数微分方程系数回归
您好斑竹,我的最终目的是要识别k参数,我的微分方程模型为
y'=k(y-y0)*y^2+v;
(我的参数识别方法是参考一个例子,后面给贴出来。
但是这个例子是不带v输入变量的)
1我前面一段程序先假设k=3,然后利用变参量数值解法求解了(仿真)y值;(因为这个式子有个带外部输入的参数v,没法求得解析解,所以只好用数值微分方程方法求出ydata;
2我后面一段程序就是利用前边得到的t值和y值去识别参数k;
3.1生成数据xxdata为2的x数据。
functiondydt=modelsq(t,y,k,v)
dydt=k
(1)*(y-6).^k
(2)+v;
clear;
k=[-3.5483.236];
y0=6;
x=[012345791114172025303540455360708090100110120130140150]';%给出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;
fori=1:
length(t1)-1
[ty]=ode23s(@modelsq,[t1(i)t1(i+1)],y0,[],k,v(i));
y0=y(end,:
);
yy=[yy;y0];
end
xydata=[t1',yy,v'];
save('xydata.txt','xydata','-ascii');%把tt和yy存为txt文件。
3.2参数识别
方法1lsqcurvefit
functiondydt=modelsq(t,y,k,v)
dydt=k
(1)*(y-6).^k
(2)+v;
functionyi=num_value(k,x)
globalv;
y0=6;
yy=y0;
fori=1:
length(x)-1
[ty]=ode23s(@modelsq,[x(i)x(i+1)],y0,[],k,v(i));
y0=y(end,:
);
yy=[yy;y0];
end
yi=yy;
clear;
globalv;
loadxydata.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-reflective','display','iter');
lb=[-131];ub=[-15];
kp=lsqcurvefit(@num_value,k0,xdata,ydata,lb,ub,options);
%kp=[-3.5400,3.23000];
方法2multistart和lsqnonglin,全局算法(速度比较慢)
globalv;
loadxydata.txt;
xdata=xydata(:
1);
ydata=xydata(:
2);
v=xydata(:
3);
F=@(k)num_value(k,xdata)-ydata;
ms=MultiStart('Display','iter','UseParallel','always');
matlabpoolopen2;
opts=optimset('Algorithm','trust-region-reflective');
problem=createOptimProblem('lsqnonlin','x0',[-4,2],'objective',F,'lb',[-61],'ub',[-15],'options',opts);
[xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100)
matlabpoolclose
opts=optimset('Algorithm','trust-region-reflective','MaxTime',200,'Maxiter',400);
方法3lsqcurvefit和multistart
globalv;
loadxydata.txt;
xdata=xydata(:
1);
ydata=xydata(:
2);
v=xydata(:
3);
ms=MultiStart('Display','iter','UseParallel','always');
opts=optimset('Algorithm','trust-region-reflective');
matlabpoolopen
problem=createOptimProblem('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)
方法4globalsearch
globalv;
loadxydata.txt;
xdata=xydata(:
1);
ydata=xydata(:
2);
v=xydata(:
3);
F=@(k)norm(num_value(k,xdata)-ydata);
gs=GlobalSearch('Display','iter');
matlabpoolopen
opts=optimset('Algorithm','interior-point','UseParallel','always');
problem=createOptimProblem('fmincon','x0',[-4,2],'objective',F,'lb',[-61],'ub',[-15],'options',opts);
[xming,fming,flagg,outptg,manyminsg]=run(gs,problem)
方法5模拟退火
tic;
globalv;
loadxydata.txt;
xdata=xydata(:
1);
ydata=xydata(:
2);
v=xydata(:
3);
F=@(k)norm(num_value(k,xdata)-ydata);
X0=[-4,2];%Startingpoint
lb=[-61];ub=[-15];
options=saoptimset('Display','iter');
[x,fval,exitFlag,output]=simulannealbnd(F,X0,lb,ub,options)
time=toc;
二、关于求解变参数微分方程
看到有不少人问过如何处理"变参数微分方程",所以抽了一点时间,写了一个例子,
希望能帮助到那些需要求解此类问题的人。
%%%================================%%%
clearall
fun=inline('[y
(2);sin(w*t)-2*y
(1)-3*y
(2)]','t','y','flag','w');
tsp=[010];
y0=[11];
xlim(tsp)
forw=1:
10
[t,y]=ode45(fun,tsp,y0,[],w);
plot(t,y)
title(['w=',num2str(w)]);
pause(0.5);
end
三、关于求解变上限积分问题
经常看到有人询问如何求解“变上限积分问题”,现举一个例子,说明该如何求解之:
希望以后碰到类似问题时,能够自己动手解决.
%%%===============================%%%
clearall
warningoff
x=linspace(0,150);
fork=1:
length(x)
xx=x(k);
fun=inline('exp(-lam.^2)');
w(k)=0.62*2/sqrt(pi)*quadl(fun,0,sqrt(pi)*xx/150);
end
plot(x,w)
4变系数数值积分问题。
Ti=Ap*R*(p1+p3).*(sin(o)+(R*sin(2*o))./(2*L*(1-(R^2*sin(o).^2)/L^2).^(1/2)))-Ap*R*(p2+p4).*(sin(o)-(R*sin(2*o))./(2*L*(1-(R^2*sin(o).^2)/L^2).^(1/2)));%这个实际上是四缸机的扭矩
要求Ti在[0,4*pi],上的数值积分,而且每隔一固定角度间隔,要从外部输入p1,p2,p3,p4;
functiony=calTi(o,p1,p2,p3,p4)
R=52.5e-3;%曲柄半径
L=184e-3;%连杆长度
Ap=0.0071;%活塞面积
y=Ap*R*(p1+p3).*(sin(o)+(R*sin(2*o))./(2*L*(1-(R^2*sin(o).^2)/L^2).^(1/2)))-Ap*R*(p2+p4).*(sin(o)-(R*sin(2*o))./(2*L*(1-(R^2*sin(o).^2)/L^2).^(1/2)));
%matlab空间代码
clc;clear;
loadptwo.txt;
p1=ptwo(:
1);p2=ptwo(:
2);p3=ptwo(:
3);p4=ptwo(:
4);
tspan=0:
4*pi/length(p1):
4*pi;
Ti1=[];
fori=1:
length(tspan)-1
z=quadl(@calTi,tspan(i),tspan(i+1),[],[],p1(i),p2(i),p3(i),p4(i));
Ti1(i, :
)=z ;
end
5.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.forii=1:
length(kk)
5.y(ii)=f(kk(ii));
6.end
7.plot(kk,y)
复制代码
方法二(不用循环):
1.plot(linspace(0,5),arrayfun(@(k)quad(@(x) sin(k.*x).*x.^2,0,5),linspace(0,5)))
复制代码
振动方程求解
到有不少人问过二阶动力微分方程的求解问题,
现举一个简单的例子,其余的情形希望读者能举一反三,自己多思考.
%%%=========================================%%%
clearall
n=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))]',...
't','x','flag','n','M','K','F','B');
options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);
[t,x]=ode45(DL,[0,3],rand(n,1),options,n,M,K,F,B);
plot(t,x(:
1:
n))
MatLab向量化编程:
arrayfun函数的使用
自己的实验;
1t=0:
0.1:
5;T=arrayfun(@(x)x.^2,t);0-5的挨个平方
2t=0:
0.1:
5;f=@(x)x.^2;f(t);同上
3t=0:
0.1:
5;f=@(x)x.^2;y=arrayfun(f,t);同上
1求f(k)=int(sin(k*x)*x^2,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));
forii=1:
length(kk)
y2(ii)=f(kk(ii));
end
time=toc;
方法2:
向量化方法+匿名函数
tic;
y=arrayfun(@(k)quadl(@(x)sin(k*x).*x.^2,0,5),linspace(0,5));%相当于求一百组循环代入k。
time=toc;
例子1匿名函数数值微分方程
f=@(t,x)[x
(2);-t*x
(1)+exp(t)*x
(2)+3*sin(2*t)];
tspan=[3.94.0];%求解区间
y0=[28];%初值
[t,x]=ode45(f,tspan,y0);
plot(t,x(:
1),'-o',t,x(:
2),'-*');
legend('y1','y2');
title('y''''=-t*y+e^t*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多变量匿名函数
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)x^2+y^+a;
f1=f(2,3)
f1=@(a)x^2+y^+a
f2=f1(4)
f2=85
4.1匿名函数和m函数的联合
functione=myfun(x,y)
e=2.*exp(4.*x)+8.*sin(x)-y;
feval(@(x,y)myfun(x,y),2,3)
4.2使用匿名函数实现符号函数的赋值运算
>>symsx;%定义符号变量
>>z=2*x^3+4*x+5;%定义表达式
>>z1=diff(z,2)%求z的2阶导数的表达式
z1=12*x
z2=eval(['@(x)'vectorize(z1)]);%vectorize函数的功能是使内联函数适合数组运算的法则
>>z2(3)
ans=36
例子3匿名函数的求极值
f=@(x)exp(x)+x^2+x^(sqrt(x))-100;
formatlong;
x0=fzero(f,3);
x0;
f(x0);
3.1求aa在不同值时下列方程的根
f=@(a)@(x)exp(x)+x^a+x^(sqrt(x))-100;
aa=0:
0.01;2;arrayfun(@(a)fzero(f(a),4),aa);%4为搜索的初值。
例子4匿名函数在现实表示隐函数
y=@(x)fzero(@(y)(exp(y)+x^y)^(1/y)-x^2*y,1);
这样对于任意的x,只要调用y(x),就能得到相应的y值,如y
(1)
这时y只接受标量,利用arrayfun函数
Y=@(x)arrayfun(@(xxx)fz
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 参量 微分方程 处理
![提示](https://static.bdocx.com/images/bang_tan.gif)