24小时热门版块排行榜    

查看: 649  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +3 23jxep3nCNZb 2026-09-14 3/150 2026-09-15 09:38 by VMDqgnlKimtZ
[公派出国] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 6F5UbRU2I5hL 2026-09-14 3/150 2026-09-15 08:58 by C79jjtjAKjEn
[教师之家] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +3 6F5UbRU2I5hL 2026-09-14 3/150 2026-09-15 08:14 by C79jjtjAKjEn
[找工作] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +4 6F5UbRU2I5hL 2026-09-14 4/200 2026-09-15 08:02 by C79jjtjAKjEn
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +4 6F5UbRU2I5hL 2026-09-14 4/200 2026-09-15 08:02 by C79jjtjAKjEn
[考研] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 s3fFTmArrBt6 2026-09-14 4/200 2026-09-15 07:50 by C79jjtjAKjEn
[教师之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 6F5UbRU2I5hL 2026-09-14 4/200 2026-09-15 07:50 by C79jjtjAKjEn
[硕博家园] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 s3fFTmArrBt6 2026-09-14 4/200 2026-09-15 07:26 by C79jjtjAKjEn
[教师之家] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 s3fFTmArrBt6 2026-09-13 4/200 2026-09-15 06:50 by 5BDX0d0WFp7t
[硕博家园] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 LwdutQ8HoqWP 2026-09-13 4/200 2026-09-15 06:02 by 5BDX0d0WFp7t
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 3n8v2C8RimXI 2026-09-13 4/200 2026-09-15 05:50 by 5BDX0d0WFp7t
[考博] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 3n8v2C8RimXI 2026-09-13 3/150 2026-09-15 05:49 by 5BDX0d0WFp7t
[考博] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 3n8v2C8RimXI 2026-09-13 5/250 2026-09-15 05:34 by 5BDX0d0WFp7t
[考研] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +5 QUjhNVAcOSff 2026-09-13 6/300 2026-09-15 05:13 by 5BDX0d0WFp7t
[考研] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +4 QUjhNVAcOSff 2026-09-13 4/200 2026-09-15 05:00 by mibUvS8DDCwf
[找工作] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +5 3n8v2C8RimXI 2026-09-13 5/250 2026-09-14 22:13 by vZfe6xYu34yj
[博后之家] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +4 QUjhNVAcOSff 2026-09-13 4/200 2026-09-14 21:37 by vZfe6xYu34yj
[教师之家] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 LwdutQ8HoqWP 2026-09-13 4/200 2026-09-14 16:13 by Aj1rhIDL5ixY
[找工作] 浙江师范大学是怎么坑我的 +6 recruit123 2026-09-09 6/300 2026-09-14 12:00 by 萨仁其其格
[教师之家] 浙江师范大学是怎么坑我的 +10 recruit123 2026-09-09 15/750 2026-09-13 10:45 by seaskyy
信息提示
请填处理意见