| 查看: 2079 | 回复: 9 | ||
| 【悬赏金币】回答本帖问题，作者hzd250将赠送您 10 个金币 | ||
[求助]
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人回复
-
各位大佬，求一份最新Wien2k版本作为学习用途
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有157人回复
-
南科大活性流体和软物质课题组诚招2027级博士、硕士研究生和博后
已经有42人回复
-
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有24人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
hzlhm
至尊木虫 (著名写手)
- 应助: 387 (硕士)
- 金币: 17852.5
- 红花: 53
- 帖子: 2880
- 在线: 606.9小时
- 虫号: 1517335
- 注册: 2011-11-30
- 性别: GG
- 专业: 常微分方程与动力系统
hzd250
2楼2021-05-17 22:00:00
3楼2021-05-18 10:22:08
4楼2021-05-18 10:52:35
hzlhm
至尊木虫 (著名写手)
- 应助: 387 (硕士)
- 金币: 17852.5
- 红花: 53
- 帖子: 2880
- 在线: 606.9小时
- 虫号: 1517335
- 注册: 2011-11-30
- 性别: GG
- 专业: 常微分方程与动力系统
5楼2021-05-18 18:39:58
6楼2021-05-18 21:31:33
7楼2021-05-18 21:34:00
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
★ ★
独孤神宇: 金币+2, 鼓励交流 2021-05-20 21:35:31
独孤神宇: 金币+2, 鼓励交流 2021-05-20 21:35:31
|
参考下： 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 |
8楼2021-05-19 20:54:58
9楼2021-05-20 12:52:15
10楼2021-05-20 13:09:46
