| 查看: 482 | 回复: 0 | |||
[求助]
求助Matlab
|
|
请高手帮我查看一下代码有什么问题，谢谢！下面是我估算出来的初值和残差。但是残差比实验值还要大，代码有问题可以我Matlab不懂，求大神指导。 function DMM170 clear all clc k0 = [0.01 0.006 0.01 8 0.6 9 0.8 1 0.3 0.007]; lb = [0 0 0 0 0 0 0 0 0 0]; ub = [inf inf inf inf inf inf inf inf inf inf]; P0 =[0.016664 4.642001 0.340144 0 0 0; 0.019859 4.576200 0.402387 0 0 0; 0.024201 4.481438 0.492568 0 0 0; 0.017000 4.641373 0.340098 0 0 0; 0.020132 4.575695 0.402343 0 0 0; 0.024499 4.480897 0.492509 0 0 0; 0.017128 4.641134 0.340081 0 0 0; 0.020678 4.574688 0.402254 0 0 0; 0.025811 4.478524 0.492248 0 0 0; 0.046146 4.056541 0.891733 0 0 0 ]; % 初始分压,MPa Pi=[0.0052847 4.6395855 0.3479159 0.0037084 0.0032052 0.0003002; 0.0072011 4.5720394 0.4119683 0.0044129 0.0040198 0.0003585; 0.0089072 4.4774599 0.5037451 0.0048294 0.0046554 0.0004029; 0.0049652 4.6375496 0.3490244 0.0032731 0.0047272 0.0004605; 0.0063178 4.5709579 0.4129024 0.0039176 0.0053778 0.0005265; 0.0073769 4.4786796 0.5039073 0.0043511 0.0051445 0.0005405; 0.0040522 4.6367309 0.3498930 0.0030793 0.0057272 0.0005174; 0.0051847 4.5702655 0.4137574 0.0037354 0.0063233 0.0007336; 0.0062153 4.4770637 0.5052473 0.0043140 0.0064064 0.0007533; 0.0191512 4.0567809 0.9105013 0.0078227 0.0050752 0.0006686 ]; % 经过Wc/F0后，各物质分压,MPa % 使用函数lsqnonlin()进行参数估计 opt=optimset('Algorithm','levenberg-marquardt'); [k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,[],[],opt,P0,Pi); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7)) fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8)) fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9)) fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10)) fprintf('\t残差平方和 = %.4f\n',resnorm) fprintf('\texitflag = %.4f\n',exitflag) fprintf('\tresidual = %.4f\n',residual) % ------------------------------------------------------------------ function f = ObjFunc(k,P0,Pi) % 目标函数 [m,n] = size(P0); Pcal = zeros(m,n); tspan = [0 0.82; 0 0.97; 0 1.19; 0 1.02; 0 1.21; 0 1.48; 0 1.37; 0 1.62; 0 1.98; 0 2.16 ]; % 即Wc/F0,g.h/mol for i = 1:m [t PP] = ode45(@Euqations,tspan(i,  ,P0(i,,[],k);Pcal(i, = PP(end,;end f= Pcal-Pi; % ------------------------------------------------------------------ function dPdt = Euqations(t, P, k) % here t = Wc / F0 denom = 1+k(4)*P(1)+k(5)*P(3)+k(6)*P(4)+k(7)*P(5)+k(8)*P(6); % k(4) = KDMM,k(5) = KME ,k(6)=KHPM,k(7)=KPDO,k(8)=KNPA,k(9)=Kp1,k(10)=Kp2 theA =k(4)*(P(1)*P(2)-P(4)*P(3)/k(9)*P(2)) / denom^2; theB =k(6)*( P(4)*P(2)-P(5)*P(3)/k(10)*P(2))/ denom^2; theC =k(7)*P(5)*P(2)/denom^2; r1 = k(1)*theA; r2 = k(2)*theB; r3 = k(3)*theC; dPDMMdt = -r1; dPHdt = -2*r1-2*r2-r3; dPMEdt = r1+r2; dPHPMdt = r1-r2; dPPDOdt = r2-r3; dPNPAdt = r3; dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt;dPNPAdt]; 使用函数lsqnonlin()估计得到的参数值为: k1 = 0.0157 ± 5.6278 k2 = 0.0058 ± 5.8687 k3 = 0.0098 ± 462.8444 k4 = 8.0000 ± 2964.0969 k5 = 0.6002 ± 129.8991 k6 = 9.0000 ± 8644.0718 k7 = 0.8000 ± 33677.3330 k8 = 1.0000 ± 365138.2060 k9 = 0.2990 ± 15.8150 k10 = 0.0007 ± 0.2817 残差平方和 = 0.0026 exitflag = 4.0000 residual = 0.0080 residual = 0.0086 residual = 0.0103 residual = 0.0080 residual = 0.0090 residual = 0.0113 residual = 0.0082 residual = 0.0096 residual = 0.0125 residual = 0.0189 residual = -0.0044 residual = -0.0041 residual = -0.0061 residual = -0.0042 residual = -0.0049 residual = -0.0094 residual = -0.0054 residual = -0.0073 residual = -0.0127 residual = -0.0164 residual = -0.0044 residual = -0.0055 residual = -0.0062 residual = -0.0049 residual = -0.0057 residual = -0.0056 residual = -0.0049 residual = -0.0056 residual = -0.0059 residual = -0.0107 residual = -0.0003 residual = -0.0003 residual = 0.0002 residual = 0.0007 residual = 0.0009 residual = 0.0014 residual = 0.0018 residual = 0.0021 residual = 0.0027 residual = 0.0002 residual = -0.0032 residual = -0.0040 residual = -0.0046 residual = -0.0047 residual = -0.0054 residual = -0.0051 residual = -0.0057 residual = -0.0063 residual = -0.0064 residual = -0.0051 residual = -0.0003 residual = -0.0004 residual = -0.0004 residual = -0.0005 residual = -0.0005 residual = -0.0005 residual = -0.0005 residual = -0.0007 residual = -0.0008 residual = -0.0007 >> |
回复此楼» 猜你喜欢
【复旦大学】二维材料方向招收2026年博士研究生1名
已经有0人回复
-
北京纳米能源与系统研究所 王中林院士/曹南颖研究员课题组2026级硕/博/博后招生
已经有10人回复
-
物理学I论文润色/翻译怎么收费?
已经有265人回复
-
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生
已经有23人回复
-
反铁磁体中的磁性切换：两种不同的机制已成功可视化
已经有0人回复
-
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
-
新西兰Robinson研究所招收全奖PhD
已经有0人回复
-
石墨烯转移--二氧化硅衬底石墨烯
已经有0人回复