24小时热门版块排行榜    

查看: 638  |  回复: 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的回帖

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的回帖

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的回帖
相关版块跳转 我要订阅楼主 laojin 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 8/400 2026-08-29 08:32 by foolishmani
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +13 晴天加油 2026-08-26 14/700 2026-08-28 18:36 by huagongfeihu
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +4 gdd2018 2026-08-28 9/450 2026-08-28 14:22 by jklily
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 申请删除本帖 +6 lyz123lyz 2026-08-27 7/350 2026-08-27 17:31 by 宁静致远sy
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 怎么看青基中了没有啊 +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
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +7 process2012 2026-08-23 10/500 2026-08-26 19:23 by hmhminy
[基金申请] 范进中举一文的中心思想 +9 炎黄贵胄 2026-08-22 10/500 2026-08-26 15:40 by semaglutide
[基金申请] 能否退出参与的面上项目解除限项 +23 koalala 2026-08-24 26/1300 2026-08-26 14:29 by 宝贝虫子
[基金申请] 出来了 +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
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+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 南海小哥
信息提示
请填处理意见