24小时热门版块排行榜    

查看: 639  |  回复: 2
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

laojin

木虫 (小有名气)

[求助] 用matlab求动力学参数代码有错误 已有1人参与

初学者,所以编了一个最简单的方程,但是有一点问题。
反应方程:-dc/dt=k*c^a

function KA
format long
clear all
clc
tspan = [0  2  4  6 8 10];
x0 = 0.846812;
k0 = [0.02  0];   
lb = [0  -1];
ub = [1  1];

data=...
    [

0        0.846812
2        0.802778
4        0.757897
6        0.702854
8        0.653739
10        0.609704

];
yexp = data(:,2);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)

function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
ysim(:,1) = Xsim1(1:end);

size(ysim(:,1));
size(yexp(:,1));

f = (ysim(:,1)-yexp(:,1));


function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
dCAdt = -k(1)*(C^k(2));
dCdt = dCAdt;

运行之后,得到结果
使用函数lsqnonlin()估计得到的参数值为:
        k1 = 0.021862398 ± 0.005117281
        k2 = -0.282935990 ± 0.797237745
  The sum of the squares is: 4.866526201e-05
K1还好,但是K2显然有问题,我是按照论坛里的Code改的,不知道出了什么问题,请各位帮忙看看。
我自己编了一个零级反应的数据,但是也是K2出现问题。
不知道是不是数据太少的缘故。
烦请大神修正,不胜感激。
回复此楼

» 猜你喜欢

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

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

laojin

木虫 (小有名气)

引用回帖:
2楼: Originally posted by dingd at 2015-01-23 10:34:55
应该没错吧,1stOpt求解:

Variable t,c;
ODEFunction c'=-k*c^a;
Data;
0        0.846812
2        0.802778
4        0.757897
6        0.702854
8        0.653739
10        0.609704


均方差 ...

只是觉得K2的偏差太大了,不知道是不是数据量太少的缘故。
3楼2015-01-23 11:44:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 3 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
laojin: 金币+10, ★★★很有帮助 2015-01-23 11:44:11
应该没错吧,1stOpt求解:
CODE:
Variable t,c;
ODEFunction c'=-k*c^a;
Data;
0        0.846812
2        0.802778
4        0.757897
6        0.702854
8        0.653739
10        0.609704

均方差(RMSE):0.00311978140170499
残差平方和(SSE):4.86651799721219E-5
相关系数(R): 0.999116092954796
相关系数之平方(R^2): 0.998232967201257
确定系数(DC): 0.997978305044057
F统计(F-Statistic): 1483.90339065772

参数                  最佳估算
--------------------        -------------
k        0.0218659845981955
a        -0.282434977032699

====== 结果输出 ======

文件: 数据文件-1
No        t        目标 c        计算 c
1        2        0.802778        0.800616656747102
2        4        0.757897        0.753655380090554
3        6        0.702854        0.705852003645687
4        8        0.653739        0.657115339412452
5        10        0.609704        0.607334455281825
2楼2015-01-23 10:34:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 9/450 2026-08-29 09:59 by jklily
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 09:39 by G6APbkg8SA6w
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 ASdOkHsho7FD 2026-08-28 7/350 2026-08-29 09:38 by G6APbkg8SA6w
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +12 Kittylucky 2026-08-27 13/650 2026-08-29 00:04 by superceng
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +4 gdd2018 2026-08-28 9/450 2026-08-28 14:22 by jklily
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 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 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见