当前位置: 首页 > 计算模拟 >matlab拟合反应动力学

matlab拟合反应动力学

作者 hzd250
来源: 小木虫 450 9 举报帖子
+关注

刚入门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比 返回小木虫查看更多

今日热帖
  • 精华评论
  • hzlhm

    有具体的数据吗?可以发给我吗?,可以试一试帮你找问题?

  • hzlhm

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

    ...

    每列数据对应的变量是什么?

  • hzd250

    每列数据对应的变量就是五中物质的浓度变化:x1 x2 x3 x4 x5

  • dingd

    引用回帖:
    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

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓