matlab中tolx,求助计算中matlab遇到的问题 - 计算模拟 - 小木虫 - 学术 科研 互动社区...

所得出的结果等于初始值，这个问题折磨我很久了 不知道怎么解决

求助大虾帮忙

愿意奉上BB

QQ403188374

Optimization terminated: norm of the current step is less

than OPTIONS.TolX.

global theta yexp S HI

k0 = [0.4  2.5  3.6  1000  1000  50 200 2500  5000   200  5000 10 90 ];         % 参数初值

lb = [0  0  0  0  0 0  0  0  0 0  0  0  0 ];                   % 参数下限

ub = [  +inf  +inf  +inf +inf  +inf  +inf  +inf +inf  +inf  +inf  +inf  +inf +inf ]; % 参数上限

x0=[0.000434037 0.000696735 0.001125753 0.00034319  0.000297428 0.000414565 6.96812E-06 2.79345E-07 0];

tspan=[0 1];

Kinetics =...

[1  0.00047     0.000722807 0.001132813 0.000320087 0.000290179 0.000399206 4.86E-06    9.52E-08    8.63E-06

2  0.000437    0.000721053 0.001149219 0.000331224 0.000289286 0.000407937 4.51E-06    1.56E-07    9.24E-06

3  0.000478    0.000712281 0.001141406 0.000322478 0.000283036 0.000400794 5.21E-06    1.84E-07    7.41E-06

4  0.000463    0.000714912 0.001182031 0.00032137  0.000286607 0.000410635 4.98E-06    1.05E-07    8.28E-06

5  0.000483    0.000705263 0.001135625 0.00033344  0.000296429 0.000400794 3.52E-06    1.18E-07    1.17E-05

6  0.000462    0.000776316 0.001135156 0.000308688 0.000283929 0.00036746  4.03E-06    1.63E-07    1.03E-05

7  0.000491    0.000716667 0.001105469 0.000336181 0.000284821 0.000361111 3.65E-06    1.44E-07    1.13E-05

8  0.000499    0.000819298 0.001078125 0.000325277 0.000283036 0.000352381 3.36E-06    1.49E-07    1.20E-05

9  0.000495    0.000767544 0.001126563 0.000336589 0.000290179 0.000331746 4.00E-06    1.28E-07    1.05E-05

10  0.000547    0.000798246 0.001145313 0.000315831 0.000289286 0.000365873 4.00E-06    1.12E-07    1.05E-05

11  0.000509    0.000729825 0.001130313 0.000323761 0.000277679 0.000387302 3.56E-06    9.84E-08    1.17E-05

12  0.000473    0.000722807 0.001157031 0.000324257 0.000288393 0.00036619  4.91E-06    8.74E-08    8.36E-06

13  0.0005      0.000813158 0.001153125 0.000334461 0.000283036 0.000373016 1.67E-06    5.48E-08    1.64E-05

14  0.000571    0.000823684 0.001151563 0.000192187 0.000253571 0.000343651 2.77E-06    7.74E-08    1.36E-05

15  0.000546    0.000718421 0.001185938 0.000333032 0.000282143 0.000394444 3.62E-06    9.64E-08    1.14E-05

];

yexp=Kinetics(:,2:10);

theta=[259.7663376  389.6495064 519.5326752 649.415844  259.7663376 389.6495064 519.5326752 649.415844  259.7663376 389.6495064 519.5326752 649.415844  259.7663376 389.6495064 519.5326752]';

HI=[0.008232711 0.008232711 0.008232711 0.008232711 0.013572204 0.013572204 0.013572204 0.013572204 0.016198704 0.016198704 0.016198704 0.016198704 0.021367521 0.021367521 0.021367521]';

S=[1.525111607  2.033482143 2.541852679 3.55859375  2.055803571 1.541852679 3.59765625  2.569754464 2.583705357 3.6171875   1.550223214 2.066964286 3.65625 2.611607143 2.089285714]';

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@ZObjFunc,k0,lb,ub,[],tspan,x0,theta,S,HI,yexp)

function f = ZObjFunc(k,tspan,x0,theta,S,HI,yexp)           % 目标函数

for i=1:length(theta)

[t Xsim] = ode45(@Zhydrosulf,tspan,x0,[],k,theta(i),S(i),HI(i));

ysim(i,1) = Xsim(end,1);

ysim(i,2) = Xsim(end,2);

ysim(i,3) = Xsim(end,3);

ysim(i,4) = Xsim(end,4);

ysim(i,5) = Xsim(end,5);

ysim(i,6) = Xsim(end,6);

ysim(i,7) = Xsim(end,7);

ysim(i,8) = Xsim(end,8);

ysim(i,9) = Xsim(end,9);

end

f = [ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2);  ysim(:,3)-yexp(:,3);ysim(:,4)-yexp(:,4); ysim(:,5)-yexp(:,5);  ysim(:,6)-yexp(:,6);ysim(:,7)-yexp(:,7); ysim(:,8)-yexp(:,8);  ysim(:,9)-yexp(:,9)];

function dCdt = Zhydrosulf (t,C,k,theta,S,HI)              % ODE模型方程

dens = 1/(C(1)+C(2)+C(3)+C(4)+C(5)+C(6)+0.000641572+0.004526409+HI);

denom1 = (1+theta*dens*(k(6)*(C(4)+C(5)+C(6))+k(7)*0.004526409+k(8)*(C(7)+C(8))+k(9)*C(9)))*S;

denom2 =S*(1+theta*dens*(k(10)*(C(4)+C(5)+C(6))+k(11)*(C(7)+C(8))+(k(12)*HI)^2+k(13)*C(9)))^3

;

dCdt = ...

[theta*dens*theta*dens*k(1)*C(4)*HI / denom1

theta*dens*theta*dens*k(2)*C(5)*HI / denom1

theta*dens*theta*dens*k(3)*C(6)*HI / denom1

-theta*dens*theta*dens*k(1)*C(4)*HI / denom1

-theta*dens*theta*dens*k(2)*C(5)*HI / denom1

-theta*dens*theta*dens*k(3)*C(6)*HI / denom1

-theta*dens*theta*dens*k(4)*C(7)*HI / denom2

-theta*dens*theta*dens*k(5)*C(8)*HI / denom2

theta*dens*theta*dens*k(4)*C(7)*HI / denom2+theta*dens*theta*dens*k(5)*C(8)*HI / denom2];

[Last edited by billduke on 2011-4-10 at 17:11]