24小时热门版块排行榜

返回列表

wodaifei

银虫 (小有名气)

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

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

各位师哥师姐，小弟求助：
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)
可是运行不出来，小弟是新手，不知道错了哪里，希望各位高手帮我一下，小弟真的是急死了！！！！最好能给出正确的程序，小弟刚入论坛不久，金币不多，但小弟愿拿出一半给大家一表自己的心意。谢谢了！！

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

要想拥有一切，就要让自己变得足够优秀！！！

1楼 2011-03-26 10:18:53

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

信彼南山

木虫 (著名写手)

应助: 33 (小学生)
金币: 4142.9
散金: 1221
红花: 16
帖子: 1178
在线: 233.5小时
虫号: 1133529
注册: 2010-10-27
专业: 导航、制导与传感技术

★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
xiegangmai(金币+1): 谢谢应助！ 2011-03-27 09:05:07

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);
%%%
%这个函数不对，输入参数是k，你用的是k1k2k3

赞一下(2人)

回复此楼

2楼2011-03-26 23:01:03

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wodaifei

银虫 (小有名气)

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

引用回帖:

Originally posted by 信彼南山 at 2011-03-26 23:01:03:
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);
%%%
%这个函数不对，输入参数是k，你用的是k1k2k3

那应该用什么?？？在函数里面输入k1k2k3？?

赞一下

回复此楼

要想拥有一切，就要让自己变得足够优秀！！！

3楼2011-03-28 08:19:53

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖