24小时热门版块排行榜

返回列表

wodaifei

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 3339.1
散金: 572
红花: 8
帖子: 248
在线: 351.1小时
虫号: 683543
注册: 2008-12-28
性别: GG
专业: 化工系统工程

[交流] 【求助】各位师哥师姐帮小弟检查一下复杂参数回归的程序吧！！小弟万分感谢！！！

各位师哥师姐，小弟求助：
dy1/dt=k1y1y2-k2y3
dy2/dt=k1y1y2-k2y3
dy3/dt=k2y3
y1(0)=10,y2(0)=8,y3(0)=0
小弟想用四阶-龙阁库塔积分和非线性最小二乘估计以上的动力学参数k1,k2,k3,数据如下：t       y1    y2       y3

      0       10       8       0

      20       8       7       5
      40       6       6       7
      60       5       5       8
      120    3       3       10
      180    2       2       12
      300    1       1       14
以下是小弟的程序：
k0 = [1,2,3];  %设参数初值
y0 = [10,8,0];
tspan = [0,20,40,60,120,180,300]; %反应时间
y1 = [10,8,6,5,3,2,1]';%测得的第一组数据
y2 = [8,7,6,5,3,2,1]';%测得的第二组数据
y3 = [0,5,7,8,10,12,14]';%测得的第三组数据
yexp = [y1;y2;y3];
function dy = KineticEqs(t,y,k)%动力学参数方程
dy = zeros(3,1);
dy(1)=k1*y(1)*y(2)-k2*y(3);
dy(2)=k1*y(1)*y(2)-k2*y(3);
dy(3)=k2*y(3);
function f =OptObjFunc(k,y0,yexp)  %目标函数
[t Y] = ode45(@KineticEqs,tspan,y,[],k);
f1 = y(1) - y1;
f2 = y(2) - y2;
f3 = y(3) - y3;
f = [f1; f2; f3];
[k,resnorm,resid,exitflag,output,lambda,jacobian] = ...%最小二乘法
lsqnonlin(@OptObjFunc,k0,[0 0],[],[],tspan,y0)
ci = nlparci(k,resid,jacobian)%置信区间
fprintf('Estimated Parameters:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k1)
fprintf('\tk2 = %.4f ± %.4f\n',k2)
可是运行不出来，小弟是新手，不知道错了哪里，希望各位高手帮我一下，小弟真的是急死了！！！！最好能给出正确的程序，小弟刚入论坛不久，金币不多，但小弟愿拿出一半给大家一表自己的心意。谢谢了！！

[ Last edited by wodaifei on 2011-3-26 at 10:07 ]

回复此楼