声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1815|回复: 5

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

[复制链接]
发表于 2011-1-12 18:15 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 weiniuzhu 于 2011-1-12 19:23 编辑



我最近在求解一个变参数微分方程系数识别的程序,从网上看了不少例子,本论坛关于变参数微分方程数值求解有几个例子,而关于我这方面的例子基本上没有,所以求助高手谢谢了
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
回复
分享到:

使用道具 举报

发表于 2011-1-15 16:32 | 显示全部楼层
本帖最后由 bainhome 于 2011-1-15 17:19 编辑

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.                                          Norm of      First-order
  2. Iteration  Func-count     f(x)          step          optimality   CG-iterations
  3.      0          2    2.09674e-014                          22.3
  4.      1          4    2.09674e-014    2.0181e-017           22.3            0
  5.      2          6    2.09674e-014   5.04526e-018           22.3            0

  6. Local minimum possible.

  7. lsqcurvefit stopped because the size of the current step is less than
  8. the selected value of the step size tolerance.

  9. <stopping criteria details>

  10. k1 =
  11.      3
  12. resultfx =
  13.   2.0967e-014
复制代码

点评

楼主还给我发站内短信求助,说没人理他。我这两天发烧,今天看到bainhome兄也是两个论坛追着回答了。  发表于 2011-1-15 18:28

评分

1

查看全部评分

发表于 2011-1-15 23:23 | 显示全部楼层
回复 2 # bainhome 的帖子

就是这样嘛! :handshake
我了解许多网友都是到处洒网, 以前个人会为些没学过的类别找答案, 还到处看看学习
马老弟早应帮忙出手了, 就可以让大伙清鬆些学习了
当然不必马上就回应
 楼主| 发表于 2011-1-16 11:54 | 显示全部楼层
本帖最后由 weiniuzhu 于 2011-1-16 11:57 编辑

太感谢各位了,帮我找出了原因所在,我本来还打算自己慢慢学最小二乘法呢,的确楼主的帮助帮我省了一大笔时间,特别感谢bainhome大哥,很负责人很认真的一个人(经验丰富),好像在论坛上的时间也有好多年了吧。像各位学习呵呵
 楼主| 发表于 2011-1-16 12:00 | 显示全部楼层
回复 2 # bainhome 的帖子

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

点评

楼主==发帖者!! :)  发表于 2011-1-16 14:01
 楼主| 发表于 2011-2-24 11:55 | 显示全部楼层
过年期间没法上网,回家后仔细调试了程序,找出了为什么程序结果一直得到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;
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-29 15:00 , Processed in 0.067965 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表