应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求大家帮我看看下面中matlab中的程序问题出在哪里？
111/2返回列表12下一页
查看: 1360  |  回复: 10
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

1135725495

铁杆木虫 (著名写手)

[求助] 求大家帮我看看下面中matlab中的程序问题出在哪里？ 已有1人参与

求动力学参数
方程为：r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
t          c
0        8.769713879
3        5.909461709
6        5.736047766
9        4.756936544
12        4.505799773
15        4.354160734
18        4.124025932
21        3.994725485
31        3.732949397
41        3.680613046
51        3.537445782
61        3.283577643
求k1,k2,k3,,k4
下面是写的程序，大家帮忙看看，为什么结果是直线呢？
function parafit
%  
% r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
%
% dCAdt = - r1;
clear all
clc
%        t/min   CA     / mol/L
  Kinetics=[0        7.029559717
1        3.960706476
2        3.163862085
3        2.963057756
6        2.820434159
9        2.698109429
12        2.556212898
15        2.447100967
18        2.207246045
21        2.13921916
31        2.000403189
41        1.917172016
51        1.851299168
61        1.636287388];
c=Kinetics( :,1)';
k0 = [0.671   0.02253  0.641   0.04028];         % 参数初值
lb = [0 0 0 0];                   % 参数下限
ub = [100 100  100  100];    % 参数上限
x0 = [7.029559717];
yexp = Kinetics;                  
warning off
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('ca-pr','ca-real')
x(:,1)
figure;plot(t,x);
hold on

% ------------------------------------------------------------------
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
f =  sum((y(:,2)-yexp(:,2)).^2) ;
% ------------------------------------------------------------------
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);

f1 = y(:,2) - yexp(:,2);

f = [f1];
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
Kinetics=[0        7.029559717
1        3.960706476
2        3.163862085
3        2.963057756
6        2.820434159
9        2.698109429
12        2.556212898
15        2.447100967
18        2.207246045
21        2.13921916
31        2.000403189
41        1.917172016
51        1.851299168
61        1.636287388];
c=Kinetics(:,1)';
for i=1:14
    c(i)=c(:,i);
   
r1 =(k(1)*k(2)*c(i)^2-0.25*k(3)*k(4)*(9.404-c(i))^2)/(1+k(2)*c(i)+0.5*k(4)*(9.404-c(i)));
dCAdt = - r1;
dxdt = [dCAdt];
end
结果如下：

求大家帮我看看下面中matlab中的程序问题出在哪里？
H8_DJFZ~6NZC}07J0EXL1]L.png
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-05-25 21:08:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
1135725495: 金币+10, ★有帮助 2015-05-26 15:35:53
两组参考结果：

Correlation Coef. (R): 0.986609345793065
R-Square: 0.97339800120622
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093829
Parameters Best Estimate
-------------------- -------------
k1 0.0633244071641639
k2 -0.123667552854961
k3 0.0140433904054541
k4 -0.681989563148558
====== Output Results =====
File: Data file - 1
No Obs. c Cal. c
1 5.909461709 6.03887691479301
2 5.736047766 5.36051185683508
3 4.756936544 4.91507500618295
4 4.505799773 4.59634816701972
5 4.354160734 4.35722735701414
6 4.124025932 4.17236164824689
7 3.994725485 4.02646466838814
8 3.732949397 3.71546906582311
9 3.680613046 3.55315135651946
10 3.537445782 3.46435427377623
11 3.283577643 3.41460522937653

Correlation Coef. (R): 0.98660934579373
R-Square: 0.973398001207533
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093831
Parameters Best Estimate
-------------------- -------------
k1 0.481614078106072
k2 1.61663483596038
k3 0.0205007668163867
k4 46.4476555663092
====== Output Results =====
File: Data file - 1
No Obs. c Cal. c
1 5.909461709 6.03887691237705
2 5.736047766 5.36051185622609
3 4.756936544 4.91507500636805
4 4.505799773 4.5963481675709
5 4.354160734 4.35722735771124
6 4.124025932 4.17236164896553
7 3.994725485 4.02646466905497
8 3.732949397 3.71546906608907
9 3.680613046 3.55315135632873
10 3.537445782 3.46435427319561
11 3.283577643 3.41460522849613
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
2楼2015-05-26 13:56:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)
引用回帖:
2楼: Originally posted by 月只蓝 at 2015-05-26 13:56:42
两组参考结果：

Correlation Coef. (R): 0.986609345793065
R-Square: 0.97339800120622
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093829
Parameters Best Estimate
--------- ...

您好，您的计算结果我看到了，我想说的是我上面编的Matlab程序正确吗？为什么会出现那样的计算结果？还有就是您计算的方法是最小二乘吗？我现在需要用matlab计算，因为写论文需要。所以您能不能帮忙看一下程序是否正确？
赞一下 回复此楼
3楼2015-05-26 15:38:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
引用回帖:
3楼: Originally posted by 1135725495 at 2015-05-26 15:38:21
您好，您的计算结果我看到了，我想说的是我上面编的Matlab程序正确吗？为什么会出现那样的计算结果？还有就是您计算的方法是最小二乘吗？我现在需要用matlab计算，因为写论文需要。所以您能不能帮忙看一下程 ...

拟合曲线和散点差异明显，可能是因为参数初值选取得不好。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
4楼2015-05-26 16:25:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)
引用回帖:
3楼: Originally posted by 1135725495 at 2015-05-26 15:38:21
您好，您的计算结果我看到了，我想说的是我上面编的Matlab程序正确吗？为什么会出现那样的计算结果？还有就是您计算的方法是最小二乘吗？我现在需要用matlab计算，因为写论文需要。所以您能不能帮忙看一下程 ...

程序有错，最后面函数写错了

function parafit
%  
% r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
%
% dCAdt = - r1;
clear all
clc
%        t/min   CA     / mol/L
  Kinetics=[0        7.029559717
1        3.960706476
2        3.163862085
3        2.963057756
6        2.820434159
9        2.698109429
12        2.556212898
15        2.447100967
18        2.207246045
21        2.13921916
31        2.000403189
41        1.917172016
51        1.851299168
61        1.636287388];
k0 = [0.0671   1.253  0.02   48.028];         % 参数初值
lb = [0 0 0 0];                   % 参数下限
ub = [100 100  100  100];    % 参数上限
x0 = [7.029559717];
yexp = Kinetics;                  
warning off
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x,'o-',t,yexp(:,2),'*');legend('ca-pr','ca-real')


% ------------------------------------------------------------------
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
f =  sum((x-yexp(:,2)).^2) ;
% ------------------------------------------------------------------
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
f1 = x - yexp(:,2);
f = [f1];
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
dxdt=-(k(1)*k(2)*x^2-0.25*k(3)*k(4)*(9.404-x)^2)/(1+k(2)*x+0.5*k(4)*(9.404-x));
赞一下 回复此楼
The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
5楼2015-05-26 16:56:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)
引用回帖:
4楼: Originally posted by 月只蓝 at 2015-05-26 16:25:47
拟合曲线和散点差异明显，可能是因为参数初值选取得不好。...

改了初值了，我用别的软件计算出来的结果，用matlab计算出来的图不是一条斜线就是一条直线，感觉程序有问题。但是问题不知道出在哪里？
赞一下 回复此楼
6楼2015-05-26 17:03:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)
引用回帖:
5楼: Originally posted by dbb627 at 2015-05-26 16:56:17
程序有错，最后面函数写错了

function parafit
%  
% r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
%
% dCAdt = - r1;
clear all
clc
%        t/min   CA     / m ...

谢谢啊！请问这是改好的吗？我可以直接用吗？
赞一下 回复此楼
7楼2015-05-26 17:05:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)
引用回帖:
7楼: Originally posted by 1135725495 at 2015-05-26 17:05:15
谢谢啊！请问这是改好的吗？我可以直接用吗？...

改好，但是初值的选取不当，拟合结果不是很好。你可以试试
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

1135725495
The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
8楼2015-05-26 17:08:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)
引用回帖:
7楼: Originally posted by 1135725495 at 2015-05-26 17:05:15
谢谢啊！请问这是改好的吗？我可以直接用吗？...

还有就是我按照你写的运行了，运行不了了。
出现下面的错误
Error: File: parafit.m Line: 72 Column: 1
Function definitions are not permitted in this context.
赞一下 回复此楼
9楼2015-05-26 17:09:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)
送红花一朵
引用回帖:
8楼: Originally posted by dbb627 at 2015-05-26 17:08:44
改好，但是初值的选取不当，拟合结果不是很好。你可以试试...

恩恩，但是我把你的程序直接复制过来，怎么运行不了呢？
出现下面的问题。
Error: File: parafit.m Line: 72 Column: 1
Function definitions are not permitted in this context.
赞一下 回复此楼
10楼2015-05-26 17:13:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 1135725495 的主题更新
111/2返回列表12下一页
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 8rmuugja8q 2026-02-22 8/400 2026-02-23 12:22 by alian_214
[基金申请] 体制内长辈说体制内绝大部分一辈子在底层，如同你们一样大部分普通教师忙且收入低 +10 瞬息宇宙 2026-02-20 13/650 2026-02-23 11:23 by holypower
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 khieu8v8m0 2026-02-22 8/400 2026-02-23 09:35 by w4l55oybr1
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +5 usprnugpzw 2026-02-21 11/550 2026-02-23 09:24 by w4l55oybr1
[教师之家] 为什么中国大学工科教授们水了那么多所谓的顶会顶刊，但还是做不出宇树机器人？ +5 欢乐颂叶蓁 2026-02-21 8/400 2026-02-23 09:19 by 欢乐颂叶蓁
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 w89i99eaeh 2026-02-22 5/250 2026-02-23 08:04 by w4l55oybr1
[博后之家] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 6/300 2026-02-23 07:59 by w4l55oybr1
[博后之家] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +6 3dfhjxgsh7 2026-02-22 9/450 2026-02-23 07:49 by w4l55oybr1
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 4/200 2026-02-23 06:46 by jsjzfl
[公派出国] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 khieu8v8m0 2026-02-22 5/250 2026-02-23 06:29 by w4l55oybr1
[硕博家园] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 8/400 2026-02-23 06:24 by w4l55oybr1
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +5 3dfhjxgsh7 2026-02-22 6/300 2026-02-23 02:04 by 5jlh3qtdvx
[教师之家] 版面费该交吗 +7 苹果在哪里 2026-02-22 8/400 2026-02-22 22:37 by otani
[基金申请] 面上可以超过30页吧？ +4 阿拉贡aragon 2026-02-22 4/200 2026-02-22 21:22 by 山西悬空寺空悬�
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 usprnugpzw 2026-02-21 6/300 2026-02-22 19:48 by w89i99eaeh
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 3dfhjxgsh7 2026-02-22 4/200 2026-02-22 16:52 by khieu8v8m0
[公派出国] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 usprnugpzw 2026-02-21 4/200 2026-02-22 16:27 by khieu8v8m0
[基金申请] “人文社科而论，许多学术研究还没有达到民国时期的水平” +4 苏东坡二世 2026-02-18 5/250 2026-02-22 16:07 by liangep1573
[基金申请] 什么是人一生最重要的？ +4 瞬息宇宙 2026-02-21 4/200 2026-02-22 11:44 by huagongfeihu
[基金申请] 今年春晚有几个节目很不错，点赞！ +11 瞬息宇宙 2026-02-16 12/600 2026-02-21 21:14 by lq493392203
信息提示
关闭
请填处理意见
关闭