求助Matlab变参数常微分方程的常系数识别程序

weiniuzhu

补充资料                                

我最近在求解一个变参数微分方程系数识别的程序，从网上看了不少例子，本论坛关于变参数微分方程数值求解有几个例子，而关于我这方面的例子基本上没有，所以求助高手谢谢了

1   假设已知k=3;利用数值方法根据给出的x计算出y

首先先取k=3,利用数值方法生成变参数微分方程数值解，然后利用仿真解识别系数k。其中v为变参数，需要从外部输入。在这里我取v=1:length(x);

function dydt =modelnew(t,y,k,v)
%辅函数
global y0;
dydt =-k*(y-y0)*y.^2+v;

function y=num_valuetest(k,x)

global y0;

v=1:length(x);

yy=y0;

fori=1:length(x)-1

[t y]=ode23s(@modelnew,[x(i) x(i+1)],y0,[] ,k,v(i));

y0=y(end,: );

yy=[yy;y0];

end

y=yy;

%ttyy=[x,yy]

%save('ttyy.txt','ttyy','-ascii');%先生成数据另存为ttyy.txt
空间代码

clear all;%主程序

global y0;

y0=0;

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]'；

y=num_valuetest(3,x);%k=3;

ttyy=[x,y]

%save('ttyy.txt','ttyy','-ascii');%先生成数据另存为ttyy.txt;

2  然后利用上面仿真得到的数据，识别参数k

开始编写参数识别程序

function odetest%主程序
clc;clear;
global y0
load ttyy.txt;
xdata=ttyy(:,1);ydata=ttyy(:,2);
k0=2.7;y0 =ydata(1);
options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter');
k=lsqcurvefit(@num_value,k0,xdata,ydata,[0],[10],options);
Jc=num_value(k,xdata);
plot(xdata,ydata,'o',xdata,Jc);

functiondydt=modelsq(t,y,k,v)
global y0
dydt =-k*(y-y0)*y.^2+v;
function y=num_value(k,x)

global y0;

v=1:length(x);

yy=y0;

for i=1:length(x)-1[

t y]=ode23s(@modelsq,[x(i) x(i+1)],y0,[] ,k,v(i));
  

y0=y(end,: );

yy=[yy;y0];

end

y=yy;

得出的结果为什么k一直是初始值呢，根本没有变化，而得不到准确值3

messenger

程序有点乱呀，函数num_valuetest和函数num_value是什么关系？怎么都有微分方程求解？
save('ttyy.txt','ttyy','-ascii');怎么都注释掉了？

weiniuzhu

2# messenger
改一下我的程序（使其明了一些）：
您好斑竹，我的最终目的是要识别k参数，我的微分方程模型为
y'=k(y-y0)*y^2+v ;
(我的参数识别方法是参考一个例子，后面给贴出来。但是这个例子是不带v输入变量的）
1 我前面一段程序先假设k=3,然后利用变参量数值解法求解了（仿真）y值；(因为这个式子有个带外部输入的参数v，没法求得解析解，所以只好用数值微分方程方法求出ydata；
2 我后面一段程序就是利用前边得到的t值和y值去识别参数k；
3 两者的num-value没啥区别，但也没联系（就是前边是为仿真出y值，后边的是为参数识别lsqcurvefit准备的模型。
111111111111111   数值解出y；
前一段程序也可改为

function dydt =modelnew(t,y,v)
global y0;
dydt =-3*(y-y0)*y.^2+v;

clear all;%主程序

global y0;

y0=0;

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值

v=1:length(x);

%给出v值(外部手动输入量），在实际情况下这个v值可能是随时间变化的，但只能每隔一段时间观测一次，如果想得到精确仿真解，可以把时间细化，然后对v进行插值.

yy=y0;

fori=1:length(x)-1

[t y]=ode23s(@modelnew,[x(i) x(i+1)],y0,[],v(i));

y0=y(end,: );

yy=[yy;y0];

end

ttyy=[x,yy];

save('ttyy.txt','ttyy','-ascii');%把tt和yy存为txt文件。

222222222222 利用上面的ttyy数据识别k

functiondydt=modelsq(t,y,k,v)
global y0
dydt =-k*(y-y0)*y.^2+v;
function y=num_value(k,x)
global y0;
v=1:length(x);
yy=y0;
fori=1:length(x)-1
[t y]=ode23s(@modelsq,[x(i) x(i+1)],y0,[] ,k,v(i));
y0=y(end,: );
yy=[yy;y0];
end
y=yy;
空间代码（下面）
clc;clear;
global y0
load ttyy.txt;
xdata=ttyy(:,1);
ydata=ttyy(:,2);
k0=2.7;
y0 =ydata(1);
options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter');
k=lsqcurvefit(@num_value,k0,xdata,ydata,[0],[10],options);%这个lsqcurvefit的主要原理是调整参数k（根据k0初始值)使得ydata与数值求解出的（num_value函数得出的）残差不断减小。

weiniuzhu

2# messenger
这是我想出的一种近似解法，就是对y拟合然后进行差分求出y'数据，把微分方程变为代数方程
y'=-k(y-y0)*y^2+v ；这是我要识别微分方程
（令dy/dt=z,;y=x1;v=x2;这样就变为一个多变量参数识别模型：z=-k*x1^3+x2;）
这里我用另一种方法去识别这个变参数微分方程的系数（但我认为这不是一种特别理想的方法，所以如果另有高手对这个问题有新的思路，请指教，我的qq393773310)，就是对仿真出的y值进行数值差分,diff方法求得y',这样上述方程就变为一个代数方程，这样就可以直接用lsqcurvefit方法求参数。
为了简单起见，我把上边的模型改为y'=-ky^3+v;
第一步：首先利用数值方法求解y数据，为下面参数识别做准备
function dydt =modelnew(t,y,v)
dydt =-3*y.^3+v;
空间代码
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]';
t1=min(x):(max(x)-min(x))/1000:max(x);
v=1:length(x);
v1=spline(x,v,t1) ;
y0=0;
yy=y0;
for i=1:length(t1)-1
[t y]=ode23s(@modelnew,[t1(i) t1(i+1)],y0,[],v1(i));
y0=y(end,:);
yy=[yy;y0];
end

z1=diff(yy)/ ((max(x)-min(x))/1000);%对yy数据进行差分，并除以Δt
y1=yy(1:length(t1)-1);%由于差分后，数据个数减少1个。
v1= v1(1:length(t1)-1);
zxxdata=[y1,v1',z1];
save('testdata.txt','zxxdata','-ascii');%存为txt文件
第二步：利用上面数据进行参数识别
%令dy/dt=z,;y=x1;v=x2;这样就变为一个多变量参数识别模型：z=-k*x1^3+x2;
2.1 模型M函数

functiony=dlsq(k,xdata)
y=-k*xdata(:,1).^3+xdata(:,2);
2.2 %空间代码
load testdata.txt;
xdata=testdata(:,1:2);ydata=testdata(:,3);
k0=[1.8]';
options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',10000,'Algorithm','trust-region-reflective');
beta=lsqcurvefit(@dlsq,k0,xdata,ydata,[],[],options);
结果k=3.0010;
这种方法看起来是有效地，不过实际还有待检验，可能对数据y进行差分的过程，如果时间间隔取得太大，那么y'数据将不准，会对参数识别k影响较大，所以只有加大数据量，这个可能是弊端吧

weiniuzhu

2# messenger
这是我参考的微分方程参数识别例子
问题如下：
已知实验数据x,y，并且x,y的关系满足以下常微分方程
Dy/dx=-k*(y-y0)*y^2
其中 k是需要回归的参数，y0是一个常数，通常等于y向量中的最后一个数值。要求：
1.通过lsqcurvefit或者lsqnonlin回归出系数k
2.画出模型预测值和实验值的对比图，模型预测值可以通过得到常微分方程数值解后三次样条spline插值得到。我已经写好的程序如下，里面有错误，我自己找不出来，请高手帮帮忙，谢谢啊
可以加我的QQ交流：40231185
=======================================
function odetest
clc;clear;
global Je J0
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);
function y=cakefun(beta,x)
global J0
tspan=[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 Je
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

bainhome

刚才发了个帖子，又删掉了，因为你四楼给的程序x赋值出现错误，后面少了个转置，改完程序能跑。
但是原建议不变，这个global用得极烂，多个function中出现clear和对y0的数值、维数动态改变，感觉这种参数传递是导致出现结果不变的首要症结。

weiniuzhu

谢谢大虾帮忙，我再试试，

bainhome

收回刚才的建议，还是没仔细跑，粗看还真是不好使，抱歉。
你的问题中主要错误：

1. options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter');

复制代码

这个optimset的设置纯属扯淡，建议如下：
1.双精度格式类型变量最高精度只能到1e-16，设到1e-20无任何意义
2.函数容差TolFun、TolX的设定是个技术活儿，一般容差设定在-4到-5的精度已经很不错了，lb和ub的选择范围太大，如果0到10之间都可以，那么lsqcurvefit会认为你给出的要求不严，很快达收敛条件
3.函数最大执行次数改为1000
4.算法自己不用设定。
odetest更改代码：

1. function odetest
   
2. %主程序
   
3. %%
   
4. clc;
   
5. global y0
   
6. load ttyy.txt;
   
7. xdata=ttyy(:,1);
   
8. ydata=ttyy(:,2);
   
9. k0=400;
   
10. y0 =ydata(1);
    
11. options=optimset('TolFun',1e-6,'TolX',1e-4,'MaxFunEvals',1000,'Display','iter');
    
12. [k1,resultfx]=lsqcurvefit(@num_value,k0,xdata,ydata,2,4,options);
    
13. Jc=num_value(k1,xdata);
    
14. plot(xdata,ydata,'o',xdata,Jc);

复制代码

注意初始k0取很大都没关系，但是TolFun、TolX、lb和ub的取值就比较精妙，需要自己更改调整。
我的用

1. k0=4000
   
2. [lb,ub]=[2,5]

复制代码

条件得到结果：

1. 
   
2.                                          Norm of      First-order
   
3. Iteration  Func-count     f(x)          step          optimality   CG-iterations
   
4.      0          2         1.51425                     4.73e+008
   
5.      1          4         1.51425    1.8248e-009      4.73e+008            0
   
6.      2          6         1.51425   4.56199e-010      4.73e+008            0
   
7. 
   
8. Local minimum possible.
   
9. 
   
10. lsqcurvefit stopped because the size of the current step is less than
    
11. the selected value of the step size tolerance.
    
12. 
    
13. <stopping criteria details>
    
14. 
    
15. k1 =
    
16.     3.5000
    
17. resultfx =
    
18.     1.5142

复制代码

用

1. k0=4000
   
2. [lb,ub]=[1,4]

复制代码

条件得到结果：

1. 
   
2.                                          Norm of      First-order
   
3. Iteration  Func-count     f(x)          step          optimality   CG-iterations
   
4.      0          2          2.3676                     9.24e+008
   
5.      1          4          2.3676   1.45328e-009      9.24e+008            0
   
6.      2          6          2.3676   3.63319e-010      9.24e+008            0
   
7. 
   
8. Local minimum possible.
   
9. 
   
10. lsqcurvefit stopped because the size of the current step is less than
    
11. the selected value of the step size tolerance.
    
12. 
    
13. <stopping criteria details>
    
14. 
    
15. k1 =
    
16.     2.5000
    
17. resultfx =
    
18.     2.3676

复制代码

而用

1. k0=4000
   
2. [lb,ub]=[2,4]

复制代码

条件得到结果：

1. 
   
2.                                          Norm of      First-order
   
3. Iteration  Func-count     f(x)          step          optimality   CG-iterations
   
4.      0          2    2.09674e-014                          7.44
   
5.      1          4    2.09674e-014   3.49546e-017           7.44            0
   
6.      2          6    2.09674e-014   8.73865e-018           7.44            0
   
7. 
   
8. Local minimum possible.
   
9. 
   
10. lsqcurvefit stopped because the size of the current step is less than
    
11. the selected value of the step size tolerance.
    
12. 
    
13. <stopping criteria details>
    
14. 
    
15. k1 =
    
16.      3
    
17. resultfx =
    
18.   2.0967e-014

复制代码

以上说明MATLAB的初值依赖性很强，不想动脑子的话，用1stopt还是稍好些，但是D版无法对微分方程进行拟合，需要买一个。

bainhome

不过MATLAB也有自己的办法，对十分合理的lb和ub上下界，放一个while循环，结果也许会更妙：

1. clc;
   
2. global y0
   
3. load ttyy.txt;
   
4. xdata=ttyy(:,1);
   
5. ydata=ttyy(:,2);
   
6. k0=4000;
   
7. y0 =ydata(1);
   
8. options=optimset('TolFun',1e-4,'TolX',1e-3,'MaxFunEvals',100,'Display','iter');
   
9. k1=k0;
   
10. lb=0;ub=15;
    
11. while abs(k1-3)>=.05
    
12. [k1,resultfx]=lsqcurvefit(@num_value,k1,xdata,ydata,lb,ub,options);
    
13. if k1>3
    
14.     k1=k1-.1;
    
15. else
    
16.     k1=k1+.1;
    
17. end
    
18. lb=k1-4;
    
19. ub=k1+4;
    
20. end
    
21. k1
    

复制代码

结果：

1. 
   
2.                                          Norm of      First-order
   
3. Iteration  Func-count     f(x)          step          optimality   CG-iterations
   
4.      0          2    2.09674e-014                          22.3
   
5.      1          4    2.09674e-014    2.0181e-017           22.3            0
   
6.      2          6    2.09674e-014   5.04526e-018           22.3            0
   
7. 
   
8. Local minimum possible.
   
9. 
   
10. lsqcurvefit stopped because the size of the current step is less than
    
11. the selected value of the step size tolerance.
    
12. 
    
13. <stopping criteria details>
    
14. 
    
15. k1 =
    
16.      3
    
17. resultfx =
    
18.   2.0967e-014

复制代码

weiniuzhu

9# bainhome
谢谢楼主能够在百忙中帮我调试程序，很是感激啊。我发在ilovematlab上面，没人理我，所以想起了这里，看来高手还是挺多。我也是最近才开始好好学matlab的，以前学过一些，不过没有真正做过啥课题，都学一些皮毛。看了您的程序，写的非常好，理论基础扎实，我现在就会调用matlab程序，自己编程还不行啊。看来我也得好好学学编程了，matlab 的确很强大啊，要学的东西也挺多。再次谢谢楼主了。
再麻烦一下您看matlab高效编程与应用25个案例分析中改进随机游走法能嵌套在里边吗
因为我的课题是的数据是采集卡得到的，里边含有噪声，我想让我的程序容错能力强一些，就是让参数识别更准确一些，有啥好建议吗

function[mx,minf]=reandwalk(f,x,lamda,epsilon,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:n

X(ii,:)=x+lamda*(u(ii,:)/norm(u(ii,:)));

F(ii)=f(X(ii,1),X(ii,2));

end

[f11,kk]=min(F);

if f11<f1

f1=f11;

x=X(k,:);

k=1;

else

k=k+1

end

end

lamda=lamda/2;

end

mx=X(kk,:);

minf=f1;

bainhome

这个还是去问作者吧，没理由我办事他挣钱。:lol
再说我就喜欢看他给读者回答问题，那真是百看不厌啊。

rocwoods

bainhome你也太狠啦，我都发烧了也不放过我。原计划今天晚上和大学宿舍的几个哥们好好腐败的，都没成行。楼主不着急的话等我退烧了再帮你看看吧，说实在的，现在这状态，回复的话怕烧糊涂了误导你。

bainhome

sorry，本想凑个章回本，比如什么：“qibbxxt温酒写代码，roc病中斩一城”这样的龌龊标题。
说实话这个随机行走法我还没细研究，不知道原理没发言权。不过我觉得没必要硬拿这个算法往自己课题里凑，就是凑也得分析清楚自己到底想用它来干嘛。