matlab变参量微分方程参数识别Word格式.docx
《matlab变参量微分方程参数识别Word格式.docx》由会员分享，可在线阅读，更多相关《matlab变参量微分方程参数识别Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
u1=spline（tspan,u,t1）;
v1=spline（tspan,v,t1）;
length（t1）-1
[t,y]=ode45（@fun,[t1（i）,t1（i+1）],y0,[],u1（i）,v1（i））;
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'
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;
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）;
[ty]=ode23s（@modelsq,[t1（i）t1（i+1）],y0,[],k,v（i））;
xydata=[t1'
yy,v'
];
save（'
xydata.txt'
xydata'
-ascii'
%把tt和yy存为txt文件。
3.2参数识别
方法1lsqcurvefit
functionyi=num_value（k,x）
globalv;
length（x）-1
[ty]=ode23s（@modelsq,[x（i）x（i+1）],y0,[],k,v（i））;
yi=yy;
loadxydata.txt;
xdata=xydata（:
1）;
ydata=xydata（:
2）;
v=xydata（:
3）;
k0=[4000,2000]'
%初值也可以改为k0=[-6,2]'
1e-8,'
300,'
'
display'
iter'
lb=[-131];
ub=[-15];
kp=lsqcurvefit（@num_value,k0,xdata,ydata,lb,ub,options）;
%kp=[-3.5400,3.23000];
方法2multistart和lsqnonglin,全局算法（速度比较慢）
F=@（k）num_value（k,xdata）-ydata;
ms=MultiStart（'
Display'
UseParallel'
always'
matlabpoolopen2;
opts=optimset（'
problem=createOptimProblem（'
lsqnonlin'
x0'
[-4,2],'
objective'
F,'
lb'
[-61],'
ub'
[-15],'
options'
opts）;
[xminm,fminm,flagm,outptm,manyminsm]=run（ms,problem,100）
matlabpoolclose
MaxTime'
200，'
Maxiter'
400）;
方法3lsqcurvefit和multistart
matlabpoolopen
lsqcurvefit'
[-4,2],'
@num_value,'
xdata'
xdata,'
ydata'
ydata,'
[-6,1],'
[-1,5],'
[xminm,fminm,flagm,outptm,manyminsm]=run（ms,problem,100）
方法4globalsearch
F=@（k）norm（num_value（k,xdata）-ydata）;
gs=GlobalSearch（'
interior-point'
fmincon'
[-61],'
[xming,fming,flagg,outptg,manyminsg]=run（gs,problem）
方法5模拟退火
tic;
X0=[-4,2];
%Startingpoint
lb=[-61];
options=saoptimset（'
[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'
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）;
三、关于求解变上限积分问题
经常看到有人询问如何求解“变上限积分问题”，现举一个例子，说明该如何求解之：
希望以后碰到类似问题时，能够自己动手解决．
%%%===============================%%%
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）;
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空间代码
loadptwo.txt;
p1=ptwo（:
p2=ptwo（:
p3=ptwo（:
p4=ptwo（:
4）;
4*pi/length（p1）:
4*pi;
Ti1=[];
z=quadl（@calTi,tspan（i）,tspan（i+1）,[],[],p1（i）,p2（i）,p3（i）,p4（i））;
Ti1（i,
:
）=z
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）））
复制代码
振动方程求解
到有不少人问过二阶动力微分方程的求解问题,
现举一个简单的例子,其余的情形希望读者能举一反三,自己多思考.
%%%=========================================%%%
n=3;
F=[25;
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））]'
...
'
x'
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:
T=arrayfun（@（x）x.^2,t）;
0-5的挨个平方
2t=0:
f=@（x）x.^2;
f（t）;
同上
3t=0:
f=@（x）x.^2;
y=arrayfun（f,t）;
1求f（k）=int（sin（k*x）*x^2,0,5,x）,k取[0,5]上的不同值；
方法1匿名函数
f=@（k）quadl（@（x）sin（k.*x）.*x.^2,0,5）;
%相当于多重匿名函数呵呵。
kk=linspace（0,50;
%默认分为100等分
y2=zeros（size（kk））;
forii=1:
y2（ii）=f（kk（ii））;
方法2：
向量化方法+匿名函数
tic;
y=arrayfun（@（k）quadl（@（x）sin（k*x）.*x.^2,0,5）,linspace（0,5））;
%相当于求一百组循环代入k。
例子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）;
1）,'
t,x（:
2）,'
-*'
）;
legend（'
y1'
y2'
title（'
'
=-t*y+e^t*y'
+3sin2t'
xlabel（'
ylabel（'
例子2求值
f=@（a,b,x）a*x+b;
a=1:
2;
b=2:
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;
%注意需要点（.）运算
b=10:
-1:
1;
f（a,b）
3.2匿名函数的表达式中也可以有参数的传递
b=5:
c=0.1:
0.5;
f=@（x,y）x.^2+y.^2+c;
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;
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
升级会员 