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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 理论与计算方法 » matlab拟合反应动力学
101/1返回列表
查看: 1529  |  回复: 9
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
【悬赏金币】回答本帖问题，作者hzd250将赠送您 10 个金币

hzd250

新虫 (小有名气)

[求助] matlab拟合反应动力学 已有2人参与

刚入门matlab的小白，最近想做反应的动力学，按照B站up主的视频自己写了一段代码，但是运行总是出问题：

错误使用 odearguments (第 93 行)；FUNC 必须返回列向量，出错 ode45 (第 115 行)，odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
出错 Kinetics>fun (第 67 行)，[t,x]=ode45(@func,tspan,x0,[],k)；
出错 lsqnonlin (第 218 行)，initVals.F = feval(funfcn{3},xCurrent,varargin{:});出错 Kinetics (第 18 行)，lsqnonlin(@fun,k0,lb,ub,[],yexp);%非线性最小二乘法。原因:Failure in initial objective function evaluation. LSQNONLIN cannot continue.

下面是我写的代码，读取的Excel表格里有5列*7行的实验数据，劳烦大佬帮我瞅瞅哪里需要改动，万分感谢。

function Kinetics
%反应一：A+B=C+M
%r=k*XA*XB-K*XC*XM
%反应二：A+C=D+M
%r=K*XA*XC-K*XD*XM
%反应三：B+C=E+M
%r=K*XB*XC-K*XE*XM
%XM=0.175
clc
clear all;
global a b
tspan=[0.5 1 4 6 8 12 16];
yexp=xlsread('reaction.xls');
k0=[0.1 0.01 0.01 0.001 0.001 0.001];%参数初值
lb=[0 0 0 0 0 0];%下边界
ub=[+inf +inf +inf +inf +inf +inf];%上边界
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
    lsqnonlin(@fun,k0,lb,ub,[],yexp);%非线性最小二乘法
tspan=[0.5 1 4 6 8 12 16];
a=1;
b=a+6;
x0=yexp(a,;%积分初值
[t,x]=ode45(@func,tspan,x0,[],k);
t1=linspace(0.5,16,200);
ya1=spline(t,x(:,1),t1);%动力学计算得到的点进行样条插值
ya2=spline(t,x(:,2),t1);
ya3=spline(t,x(:,3),t1);
ya4=spline(t,x(:,4),t1);
ya5=spline(t,x(:,5),t1);
for m=1:7
    for n=1:5
        yy(a+m-1,n)=x(m,n);%每一次的值存入yy矩阵
    end
end
figure(1)
plot(tspan,yexp(a:b,1),'k^',t1,ya1,'k-',tspan,yexp(a:b,2),'ro',t1,ya2,'r-',tspan,yexp(a:b,3),'bd',t1,ya3,'b-',...
tspan,yexp(a:b,4),'g*',t1,ya4,'g-',tspan,yexp(a:b,5),'yp',t1,ya5,'y-');
legend('','A浓度','','B浓度','','C浓度','','D浓度','','E浓度');
xlabel('t(h)');ylabel('浓度(mol/L)');title('170℃ 0.1wt%催化剂');
t1=linspace(0.5,16,200);
z1=spline(t,yy(1:7,1),t1);
h1=spline(t,yy(1:7,2),t1);
s1=spline(t,yy(1:7,3),t1);
b1=spline(t,yy(1:7,4),t1);
u1=spline(t,yy(1:7,5),t1);
xlswrite('result.xls',[t1' z1' h1' s1' b1' u1'],'sheet1');
xlswrite('result.xls',residual,'sheet2');
Ne = length(yexp(:,2));     %模型适定性判别
Np = length(k);
[rho2,F] = rho2_F(k,yexp,resnorm,Ne,Np);
ci=nlparci(k,residual,jacobian)
fprintf('\t k1,0=%.1f ± %.4f\n',k(1),ci(1,2)-k(1));
fprintf('\t k2,0=%.1f ± %.4f\n',k(2),ci(2,2)-k(2));
fprintf('\t k3,0=%.1f ± %.4f\n',k(3),ci(3,2)-k(3));
fprintf('\t k4,0=%.1f ± %.4f\n',k(4),ci(4,2)-k(4));
fprintf('\t k5,0=%.1f ± %.4f\n',k(5),ci(5,2)-k(5));
fprintf('\t 残差平方和：%.3f\n',resnorm)
fprintf('\t 实验点数和自由度分别为 Ne = %d和 Np = %d\n',Ne,Np)
fprintf('\t 决定性指标ρ^2: %.4f\n',rho2)
fprintf('\t F比: %.3f\n\n',F)
%=================================================================================
function f=fun(k,yexp)
f=[];
tspan=[0.5 1 4 6 8 12 16];
a=1;
x0=yexp(a,;
[t,x]=ode45(@func,tspan,x0,[],k)
d=a+6;
yc1=x(:,1);
yc2=x(:,2);
yc3=x(:,3);
yc4=x(:,4);
yc5=x(:,5);
f11=yexp(a:d,1)-yc1;
f12=yexp(a:d,2)-yc2;
f13=yexp(a:d,3)-yc3;
f14=yexp(a:d,4)-yc4;
f15=yexp(a:d,5)-yc5;
ff=[f11 f12 f13 f14 f15];
f=[f;ff];
%=================================================================================
function dxdt=func(t,x,k)
r1=-k(1)*x(1)*x(2)-k(2)*x(1)*x(3)+k(4)*x(3)*0.175+k(5)*x(4)*0.175;
r2=-k(1)*x(1)*x(2)-k(3)*x(2)*x(3)+k(4)*x(3)*0.175+k(6)*x(5)*0.175;
r3=k(1)*x(1)*x(2)+k(5)*x(4)*0.175+k(6)*x(5)*0.175-k(2)*x(1)*x(3)-k(3)*x(2)*x(3)-k(4)*x(3)*0.175;
r4=k(2)*x(1)*x(3)-k(5)*x(4)*0.175;
r5=k(3)*x(2)*x(3)-k(6)*x(5)*0.175;
dxdt=[r1 r2 r3 r4 r5]
%=================================================================================
function [rho2,F] = rho2_F(k,yexp,s,Ne,Np)
y=yexp.^2;
sy = sum(y();
rho2 = 1 - s/sy;              %rho2: 决定性指标
F = (sy - s)*(Ne-Np)/(Np*s);  %F：F比
回复此楼

» 猜你喜欢
1楼 2021-05-12 19:59:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

有具体的数据吗？可以发给我吗？，可以试一试帮你找问题？
赞一下(1人) 回复此楼

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

hzd250
QQ：2120156492
2楼2021-05-17 22:00:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
送红花一朵
引用回帖:
2楼: Originally posted by hzlhm at 2021-05-17 22:00:00
有具体的数据吗？可以发给我吗？，可以试一试帮你找问题？

有的，太谢谢你了
matlab拟合反应动力学



发自小木虫Android客户端
赞一下 回复此楼
3楼2021-05-18 10:22:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
引用回帖:
3楼: Originally posted by hzd250 at 2021-05-18 10:22:08
有的，太谢谢你了

...

这个是简化的数据
matlab拟合反应动力学-1



发自小木虫Android客户端
赞一下 回复此楼
4楼2021-05-18 10:52:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

引用回帖:
4楼: Originally posted by hzd250 at 2021-05-18 10:52:35
这个是简化的数据

...

每列数据对应的变量是什么？
赞一下(1人) 回复此楼

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

hzd250
QQ：2120156492
5楼2021-05-18 18:39:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
每列数据对应的变量就是五中物质的浓度变化：x1 x2 x3 x4 x5

发自小木虫Android客户端
赞一下 回复此楼
6楼2021-05-18 21:31:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
送红花一朵
引用回帖:
5楼: Originally posted by hzlhm at 2021-05-18 18:39:58
每列数据对应的变量是什么？...

每列数据对应的变量就是五中物质的浓度变化：x1 x2 x3 x4 x5，就是后面的微分速率方程里面的x1 x2 x3 x4 x5

发自小木虫Android客户端
赞一下 回复此楼
7楼2021-05-18 21:34:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★
独孤神宇: 金币+2, 鼓励交流 2021-05-20 21:35:31
引用回帖:
4楼: Originally posted by hzd250 at 2021-05-18 10:52:35
这个是简化的数据

...

参考下：

Root of Mean Square Error (RMSE): 0.0602102835008095
Sum of Squared Residual: 0.108758347177436
Correlation Coef. (R): 0.963380309380527
R-Square: 0.928101620502121

Parameter                  Best Estimate
--------------------        -------------
k1        0.0431976858691231
k2        0.0204925345429154
k4        8.27595775884273E-20
k5        3.49573492337918E-16
k3        0.0313036416446144
k6        2.17399315624263E-18
赞一下 回复此楼

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

hzd250 hzd250
8楼2021-05-19 20:54:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
送红花一朵
引用回帖:
8楼: Originally posted by dingd at 2021-05-19 20:54:58
参考下：

Root of Mean Square Error (RMSE): 0.0602102835008095
Sum of Squared Residual: 0.108758347177436
Correlation Coef. (R): 0.963380309380527
R-Square: 0.928101620502121

Parameter       ...

太谢谢啦！
大佬，输出结果怎么才能用科学记数法表示呢，我这边输出的数据都是小数表示的

发自小木虫Android客户端
赞一下 回复此楼
9楼2021-05-20 12:52:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzd250

新虫 (小有名气)
送红花一朵
引用回帖:
8楼: Originally posted by dingd at 2021-05-19 20:54:58
参考下：

Root of Mean Square Error (RMSE): 0.0602102835008095
Sum of Squared Residual: 0.108758347177436
Correlation Coef. (R): 0.963380309380527
R-Square: 0.928101620502121

Parameter       ...

大佬，我运行的k1 k2 k3值都和你给的一样，但是后面三个k值都是2.22×10-4这是因为什么原因呢

发自小木虫Android客户端
赞一下 回复此楼
10楼2021-05-20 13:09:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 hzd250 的主题更新
101/1返回列表
不应助 确定回帖应助 (注意：应助才可能被奖励，但不允许灌水，必须填写15个字符以上)
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭