求助高手帮忙调试变参数常微分方程的常系数识别程序

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 12 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

bainhome

simwe上回答过，不过响应chaqing老哥“在振动现身”的光荣号召，过来把问题再说一遍：
问题出在：

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

复制代码

所给初值k0=2.7本身就在[0,10]之间，初值就是收敛的，当然最终结果就好像没有运算一样。
另外，很多人喜欢把TolFun和TolX值设定得很小，这位同学居然设定到1e-20，这是一种非常扯蛋的做法，原因很简单：double类型才1e-16的容差，设到20计算机能认识？
解决方案：
初值可以加到很大，但是lu和ub要缩小范围，可以把lb和ub与计算的k值结合起来，不断通过while进行调整，最终会得到比较满意的效果，这是最终进行搜索的关键之一，另外TolFun这些就不用说了吧。

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

复制代码

当然，加个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
    

复制代码

ChaChing

回复 2 # bainhome 的帖子

就是这样嘛! :handshake
我了解许多网友都是到处洒网, 以前个人会为些没学过的类别找答案, 还到处看看学习
马老弟早应帮忙出手了, 就可以让大伙清鬆些学习了
当然不必马上就回应

weiniuzhu

太感谢各位了，帮我找出了原因所在，我本来还打算自己慢慢学最小二乘法呢，的确楼主的帮助帮我省了一大笔时间，特别感谢bainhome大哥，很负责人很认真的一个人（经验丰富），好像在论坛上的时间也有好多年了吧。像各位学习呵呵

weiniuzhu

回复 2 # bainhome 的帖子

非常感谢楼主，向楼主学习

weiniuzhu

过年期间没法上网，回家后仔细调试了程序，找出了为什么程序结果一直得到3的原因，果然如同bainhome老兄在simwe论坛上所言，问题出在function y=num_value(k,x)这个M函数中的global  y0这句上，看来global还真是不能乱用，下面是修改后的程序

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文件。

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

空间代码

clc;clear;
load ttyy.txt;
xdata=ttyy(:,1);
ydata=ttyy(:,2);
k0=2.7;
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;修改前，参数传递可能有问题

y0=0;%修改后，直接改为明确的数，或者使用y=num_value(k,x,y0)来传递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;