求助高手帮忙调试变参数常微分方程的常系数识别程序
本帖最后由 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=ode23s(@modelnew,,y0,[] ,k,v(i)); y0=y(end,: );yy=;endy=yy;%ttyy=%save('ttyy.txt','ttyy','-ascii');%先生成数据另存为ttyy.txt
空间代码clear all;%主程序global y0;y0=0;x=；y=num_valuetest(3,x);%k=3;ttyy=%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,,,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,,y0,[] ,k,v(i));
y0=y(end,: );
yy=;
end
y=yy; 得出的结果为什么k一直是初始值呢，根本没有变化，而得不到准确值3 本帖最后由 bainhome 于 2011-1-15 17:19 编辑
simwe上回答过，不过响应chaqing老哥“在振动现身”的光荣号召，过来把问题再说一遍：
问题出在：options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter');所给初值k0=2.7本身就在之间，初值就是收敛的，当然最终结果就好像没有运算一样。
另外，很多人喜欢把TolFun和TolX值设定得很小，这位同学居然设定到1e-20，这是一种非常扯蛋的做法，原因很简单：double类型才1e-16的容差，设到20计算机能认识？
解决方案：
初值可以加到很大，但是lu和ub要缩小范围，可以把lb和ub与计算的k值结合起来，不断通过while进行调整，最终会得到比较满意的效果，这是最终进行搜索的关键之一，另外TolFun这些就不用说了吧。options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter');当然，加个while循环则结果更妙：clc;
global y0
load ttyy.txt;
xdata=ttyy(:,1);
ydata=ttyy(:,2);
k0=4000;
y0 =ydata(1);
options=optimset('TolFun',1e-4,'TolX',1e-3,'MaxFunEvals',100,'Display','iter');
k1=k0;
lb=0;ub=15;
while abs(k1-3)>=.05
=lsqcurvefit(@num_value,k1,xdata,ydata,lb,ub,options);
if k1>3
k1=k1-.1;
else
k1=k1+.1;
end
lb=k1-4;
ub=k1+4;
end
k1
Norm of First-order
IterationFunc-count f(x) step optimality CG-iterations
0 2 2.09674e-014 22.3
1 4 2.09674e-014 2.0181e-017 22.3 0
2 6 2.09674e-014 5.04526e-018 22.3 0
Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the selected value of the step size tolerance.
<stopping criteria details>
k1 =
3
resultfx =
2.0967e-014
回复 2 # bainhome 的帖子
就是这样嘛! :handshake
我了解许多网友都是到处洒网, 以前个人会为些没学过的类别找答案, 还到处看看学习
马老弟早应帮忙出手了, 就可以让大伙清鬆些学习了
当然不必马上就回应 本帖最后由 weiniuzhu 于 2011-1-16 11:57 编辑
太感谢各位了，帮我找出了原因所在，我本来还打算自己慢慢学最小二乘法呢，的确楼主的帮助帮我省了一大笔时间，特别感谢bainhome大哥，很负责人很认真的一个人（经验丰富），好像在论坛上的时间也有好多年了吧。像各位学习呵呵 回复 2 # bainhome 的帖子
非常感谢楼主，向楼主学习 过年期间没法上网，回家后仔细调试了程序，找出了为什么程序结果一直得到3的原因，果然如同bainhome老兄在simwe论坛上所言，问题出在function y=num_value(k,x)这个M函数中的globaly0这句上，看来global还真是不能乱用，下面是修改后的程序
clear all;%主程序global y0;y0=0;x=';%给出x值v=1:length(x);%给出v值(外部手动输入量），在实际情况下这个v值可能是随时间变化的，但只能每隔一段时间观测一次，如果想得到精确仿真解，可以把时间细化，然后对v进行插值.yy=y0;fori=1:length(x)-1=ode23s(@modelnew,,y0,[],v(i)); y0=y(end,: );yy=;endttyy=;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,,,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,,y0,[] ,k,v(i)); y0=y(end,: );
yy=;
end
y=yy;
页:
[1]