| 查看: 673 | 回复: 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 ]; |
回复此楼» 猜你喜欢
070200求调剂，一志愿某211，288分
已经有20人回复
-
法国博士后职位
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有239人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有9人回复
-
桂林理工大学物理学专业招收调剂
已经有10人回复
