| 查看: 636 | 回复: 0 | |||
[求助]
Matlab指导
|
|
这是我编的程序代码，但不知哪地方有问题，请高手帮忙指点一下。 function KineticsEstHydrogenation % clear all clc format long % KH KDMM KMe KHPM KPDO k1 k2 k3 k0 = [2.780 4.049 3.225 1.540 2.456 27.967 42.249 10]; lb = [0 0 0 0 0 0 0 0]; % 参数下限 ub = [+inf +inf +inf +inf +inf +inf +inf +inf]; % 参数上限 % 动力学数据: % T =195 % t H2in DMMin Mein HPMin PDOin NPAin HPMout PDOout NPAout data=... [2.05 0.00258 0.92932 0.06810 0 0 0 3.1360E-05 1.1642E-03 0.0011365; 1.37 0.00155 0.93029 0.06817 0 0 0 7.8438E-05 1.3477E-03 0.0005733; 1.62 0.00173 0.91758 0.08068 0 0 0 1.0444E-04 1.7271E-03 0.0004754; 1.98 0.00187 0.89928 0.09884 0 0 0 2.8942E-04 2.3631E-03 0.0004436; 1.21 0.00158 0.91772 0.08070 0 0 0 4.1373E-04 1.8968E-03 0.0002797; 1.49 0.00224 0.89895 0.09881 0 0 0 6.4058E-04 2.0193E-03 0.0002921; 2.70 0.00531 0.81544 0.17925 0 0 0 1.3040E-03 2.3203E-03 0.0003679; 1.02 0.00254 0.92936 0.06810 0 0 0 2.0611E-04 1.7761E-03 0.0003612; 0.82 0.00273 0.92918 0.06809 0 0 0 4.6464E-04 1.6084E-03 0.0001887; 0.97 0.00213 0.91722 0.08065 0 0 0 7.9002E-04 1.3764E-03 0.0001489; 1.19 0.00379 0.89756 0.09865 0 0 0 9.7162E-04 1.7297E-03 0.0002101; 2.16 0.00663 0.81436 0.17902 0 0 0 1.7438E-03 2.1190E-03 0.0003488 ]; % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[]); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1(KH) = %.4f\n',k(1)) fprintf('\tk2(KDMM) = %.4f\n',k(2)) fprintf('\tk3(KMe) = %.4f\n',k(3)) fprintf('\tk4(KHPM) = %.4f\n',k(4)) fprintf('\tk5(KPDO) = %.4f\n',k(5)) fprintf('\tk6(k1) = %.4f\n',k(6)) fprintf('\tk7(k2) = %.4f\n',k(7)) fprintf('\tk8(k3) = %.4f\n',k(8)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[]); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') Output % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[]); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') Output function f = ObjFunc4Fmincon(k) % t H2in DMMin Mein HPMin PDOin NPAin HPMout PDOout NPAout data=... [2.05 0.00258 0.92932 0.06810 0 0 0 3.1360E-05 1.1642E-03 0.0011365; 1.37 0.00155 0.93029 0.06817 0 0 0 7.8438E-05 1.3477E-03 0.0005733; 1.62 0.00173 0.91758 0.08068 0 0 0 1.0444E-04 1.7271E-03 0.0004754; 1.98 0.00187 0.89928 0.09884 0 0 0 2.8942E-04 2.3631E-03 0.0004436; 1.21 0.00158 0.91772 0.08070 0 0 0 4.1373E-04 1.8968E-03 0.0002797; 1.49 0.00224 0.89895 0.09881 0 0 0 6.4058E-04 2.0193E-03 0.0002921; 2.70 0.00531 0.81544 0.17925 0 0 0 1.3040E-03 2.3203E-03 0.0003679; 1.02 0.00254 0.92936 0.06810 0 0 0 2.0611E-04 1.7761E-03 0.0003612; 0.82 0.00273 0.92918 0.06809 0 0 0 4.6464E-04 1.6084E-03 0.0001887; 0.97 0.00213 0.91722 0.08065 0 0 0 7.9002E-04 1.3764E-03 0.0001489; 1.19 0.00379 0.89756 0.09865 0 0 0 9.7162E-04 1.7297E-03 0.0002101; 2.16 0.00663 0.81436 0.17902 0 0 0 1.7438E-03 2.1190E-03 0.0003488 ]; yexp = [data(:,8),data(:,9),data(:,10)]; for i=1:12 [t x] = ode45(@KineticEqs,[0,data(:,1)],[data(:,2),data(:,3),data(:,4),data(:,5),data(:,6),data(:,7)],[],k); y(:,1:3) = x(:,4:6) f = sum((y(:,1)-yexp(:,1)).^2+(y(:,2)-yexp(:,2)).^2+(y(:,3)-yexp(:,3)).^2); end % ------------------------------------------------------------------ function f = ObjFunc4LNL(k) % t H2in DMMin Mein HPMin PDOin NPAin HPMout PDOout NPAout data=... [2.05 0.00258 0.92932 0.06810 0 0 0 3.1360E-05 1.1642E-03 0.0011365; 1.37 0.00155 0.93029 0.06817 0 0 0 7.8438E-05 1.3477E-03 0.0005733; 1.62 0.00173 0.91758 0.08068 0 0 0 1.0444E-04 1.7271E-03 0.0004754; 1.98 0.00187 0.89928 0.09884 0 0 0 2.8942E-04 2.3631E-03 0.0004436; 1.21 0.00158 0.91772 0.08070 0 0 0 4.1373E-04 1.8968E-03 0.0002797; 1.49 0.00224 0.89895 0.09881 0 0 0 6.4058E-04 2.0193E-03 0.0002921; 2.70 0.00531 0.81544 0.17925 0 0 0 1.3040E-03 2.3203E-03 0.0003679; 1.02 0.00254 0.92936 0.06810 0 0 0 2.0611E-04 1.7761E-03 0.0003612; 0.82 0.00273 0.92918 0.06809 0 0 0 4.6464E-04 1.6084E-03 0.0001887; 0.97 0.00213 0.91722 0.08065 0 0 0 7.9002E-04 1.3764E-03 0.0001489; 1.19 0.00379 0.89756 0.09865 0 0 0 9.7162E-04 1.7297E-03 0.0002101; 2.16 0.00663 0.81436 0.17902 0 0 0 1.7438E-03 2.1190E-03 0.0003488 ]; for i=1:12 yexp = [data(:,8),data(:,9),data(:,10)]; [t x] = ode45(@KineticEqs,[0,data(:,1)],[data(:,2),data(:,3),data(:,4),data(:,5),data(:,6),data(:,7)],[],k); y(:,1:3) = x(:,4:6); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f = [f1; f2; f3]; end % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) % 动力学方程，氢气与DMO均解离吸附，表面反应为控制步骤。 r1 = k(6).*((k(1).*x(1).*k(2).*x(2)).^0.5)./... (1+2.*(k(2).*x(2)).^0.5+2.*(k(4).*x(4)).^0.5+(k(1).*x(1)).^0.5+k(3).*x(3)+k(5).*x(5)).^2; r2 = k(7).*((k(1).*x(1).*k(4).*x(4)).^0.5)./... (1+2.*(k(2).*x(2)).^0.5+2.*(k(4).*x(4)).^0.5+(k(1).*x(1)).^0.5+k(3).*x(3)+k(5).*x(5)).^2; r3 = k(8).*k(5).*k(2).*x(5).*x(2)./... (1+2.*(k(1).*x(1)).^0.5+2.*(k(4).*x(4)).^0.5+(k(2).*x(2)).^0.5+k(3).*x(3)+k(5).*x(5)).^2; dxdt=... [-2.*r1-2.*r2-r3 % H2 -r1 % DMM r1+r2 % ME r1-r2 % HPM r2-r3 % PDO r3 % NPA ]; |
回复此楼» 猜你喜欢
反铁磁体中的磁性切换：两种不同的机制已成功可视化
已经有0人回复
-
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有164人回复
-
新西兰Robinson研究所招收全奖PhD
已经有0人回复
-
石墨烯转移--二氧化硅衬底石墨烯
已经有0人回复
-
笼目材料中量子自旋液体基态的证据
已经有0人回复
-
数学教学论硕士可以读数学物理博士吗？
已经有0人回复
-
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生
已经有4人回复
-
澳门大学 应用物理及材料工程研究院 潘晖教授课题组诚招博士后
已经有11人回复