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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求大家帮我看看下面中matlab中的程序问题出在哪里？
111/2返回列表12下一页
查看: 1499  |  回复: 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下一页
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 哪位高人中了，把查询到的截图贴出来让我看看，让我长长见识 +6 yuleib84 2026-08-26 7/350 2026-08-31 19:46 by 鱼翔浅底1
[基金申请] 能否申诉？ +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[基金申请] 投票:  有多少人是今天查系统知道结果的？ +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 11:51 by xiangy672
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 国社科又开始会评了，不知道这次命运如何 +6 雨打竹帘 2026-08-30 10/500 2026-08-31 09:29 by huixian257
[基金申请] 为什么到现在没收到通知？ +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？ +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜？ 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国际合作可查了，中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 今天务委会开完了，明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 牛来！米来！面来！ +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
信息提示
关闭
请填处理意见
关闭