24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zhaoshazhu

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 125.7
散金: 10
帖子: 51
在线: 38小时
虫号: 1812564
注册: 2012-05-12
专业: 环境化工

[求助] 求Matlab高手指导已有1人参与

Matlab运算结果的置信区间很大，和什么有关系呢？
function KineticsEst6
clear all
clc
tspan = [0 662.25];
k0 = [0.4587 0.4971  10  12 6 8 205 0.5653];
lb = [0  0  0  0  0 0 0];
ub = [50  50  100  100 50 50 1000 50];

P0 =[0.04020 4.02010 0.18794 0 0;
   0.02228 4.45682 0.10418 0 0;
   0.01541 4.62428 0.07206 0 0;
   0.01336 4.67446 0.06244 0 0;
   0.04020 4.02010 0.18794 0 0;
   0.04020 4.02010 0.18794 0 0;
   0.04020 4.02010 0.18794 0 0
];  % 初始分压,MPa

Pi=[0.01303 4.00374 0.96404 0.007095 0.012097;
0.00393 4.44636 0.53626 0.074475 0.005998;
0.00127 4.61707 0.37150 0.070931 0.003074;
0.00090 4.66795 0.32206 0.066794 0.002414;
0.00243 4.00519 0.96961 0.158768 0.006883;
0.00681 4.00390 0.96766 0.117792 0.009855;
0.01303 4.00374 0.96404 0.070945 0.012097
];
% 经过Wc/F0后，各物质分压,MPa

% 使用函数lsqnonlin()进行参数估计
options=optimset('MaxFunEvals',1000000,'MaxIter',400000)
[k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,lb,ub,[],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))
% ------------------------------------------------------------------
function f = ObjFunc(k,P0,Pi)          % 目标函数
[m,n] = size(P0);
Pcal = zeros(m,n);
tspan =[0 662.25];       % 即Wc/F0,g.h/mol
for i = 1:m
[t PP] = ode45(@Euqations,tspan,P0(i,

,[],k);
Pcal(i,

= PP(end,

;
end
f= Pcal-Pi;

% ------------------------------------------------------------------
function dPdt = Euqations(t, P, k) % here t = Wc / F0
denom = 1+k(3)*P(1)+k(5)*P(4)+k(5)*P(3)+k(6)*P(5); % k(3) = KDMM, k(4) = KME ,k(5)=KHPM,k(6)=KPDO,k(7)=Kp1,k(8)=Kp2
theA =k(3)*P(1)*P(2)*(1-P(4)*P(3)/k(7)*P(1)*P(2)^2) / denom;
theB =k(5)* P(4)*P(2)*(1-P(5)*P(3)/k(8)*P(4)*P(2)^2)/ denom;
r1 = k(1)*theA;
r2 = k(2)*theB;

dPDMMdt = -r1;
dPHdt = -2*r1-2*r2;
dPMEdt = r1+r2;
dPHPMdt = r1-r2;
dPPDOdt = r2;

dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt];
使用函数lsqnonlin()估计得到的参数值为:
k1 = 0.4599 ± 89509.6910
k2 = 0.4971 ± 101889.0949
k3 = 10.0001 ± 1946900.8019
k4 = 12.0000 ± 8639183.4320
k5 = 5.9997 ± 420897.9438
k6 = 7.9999 ± 2837488.5903
k7 = 205.0000 ± 1472315356.6507
k8 = 0.5653 ± 3696268.6136
结果和我设的初值一样，是不是就没有计算程序。

回复此楼