matlab 刚性方程,MATLAB求刚性方程遇到困难请教！

有三个可逆反应组成的反应体系如下：

TG+M→←DG+ME (正逆反应速率常数分别为k1、k2)

DG+M→←MG+ME (正逆反应速率常数分别为k3、k4)

MG+M→←GL+ME (正逆反应速率常数分别为k5、k6)

建立反应以上方程的数学模型如下：

dTG/dt=-k1[TG][M]+ k2[DG][ME]

dDG/dt=k1[TG][M]- k2[DG][ME] –k3[DG][M]+ k4[MG][ME]

dMG/dt= k3[DG][M]- k4[MG][ME] –k5[MG][M]+ k6[GL][ME]

dME/dt= k1[TG][M]- k2[DG][ME]+k3[DG][M]- k4[MG][ME] +k5[MG][M]- k6[GL][ME]

其中TG、DG、MG、ME的浓度是可测定的，各物质浓度见存在如下关系：

[GL]=[TG]0-[TG]-[DG]-[MG]= 0.6380-[TG]-[DG]-[MG] ;

[M]=[M]0-[ME]=4.0626-[ME];

方程初始值为

[TG]0=0.6380

[DG]0=0

[MG]0=0

[GL]0=0

[ME]0=0

[M]0=4.0626

我的程序如下：

clear all

k0 = [120 10 30 1 5 1];

%随意给定的参数初值

lb = [0 0 0 0 0 0];

%随意给定的参数下限

ub = [1000

1000

1000

1000

1000 1000];

%随意给定的参数上限

X0=[0.638;0;0;0];

%t=0时刻四种物质TG、DG、MG、ME的初始浓度

Xexp=...

[0.6380 0 0 0

0.1512   0.1088  0.1244  1.0757

0.0955  0.0816  0.1454  1.2740

0.0719  0.0632  0.1431  1.3858

0.0500  0.0493  0.1400  1.4838

0.0451  0.0436  0.1283  1.5247

0.0361  0.0351  0.1226  1.5762

0.0346  0.0301  0.1028  1.6149

0.0261  0.0282  0.1016  1.6460

0.0237  0.0257  0.0977  1.6630]; % 25度0.8%催化剂条件下不同时间(0 0.5 1.5 2 3 4 5 7 10min)各浓度数据

%第一种计算方法，使用函数fmincon()进行参数估计

[k,fval,flag] = fmincon(@jlzengObjFmincon,k0,[],[],[],[],lb,ub,[],[],X0,Xexp);

fprintf('\n使用函数fmincon()估计得到的参数值为:\n')

fprintf('\tk1 = %.4f\n',k(1))

fprintf('\tk2 = %.4f\n',k(2))

fprintf('\tk3 = %.4f\n',k(3))

fprintf('\tk4 = %.4f\n',k(4))

fprintf('\tk5 = %.4f\n',k(5))

fprintf('\tk6 = %.4f\n',k(6))

fprintf('

The sum of the squares is: %.1e\n\n',fval)

k_fmincon = k;

%第二种计算方法，使用函数Lsqnonlin()进行参数估计

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

lsqnonlin(@jlzengObjLsqnonlin,k0,lb,ub,[],X0,Xexp);

ci=nlparci(k,residual,jacobian);

fprintf('\n\n使用lsqnonlin()函数估计得到的参数值为：\n'),Output

fprintf('\n使用函数lsqqnonlin()估计得到的参数值为:\n')

fprintf('\tk1 = %.4f\n',k(1))

fprintf('\tk2 = %.4f\n',k(2))

fprintf('\tk3 = %.4f\n',k(3))

fprintf('\tk4 = %.4f\n',k(4))

fprintf('\tk5 = %.4f\n',k(5))

fprintf('\tk6 = %.4f\n',k(6))

fprintf('

The sum of the squares is: %.1e\n\n',fval)

k_fmincon = k;

三个m文件如下：

function dXdt=jlzengkinetic(t,X,k)

%建立微分方程组，供oed23s调用

dXdt=[ (-k(1)*X(1)*(4.0626-X(4))+k(2)*X(2)*X(4))

(k(1)*X(1)*(4.0626-X(4))-k(2)*X(2)*X(4)-k(3)*X(2)*(4.0626-X(4))+k(4)*X(3)*X(4))

(k(3)*X(2)*(4.0626-X(4))-k(4)*X(3)*X(4)-k(5)*X(3)*(4.0626-X(4))+k(6)*(0.638-X(1)-X(2)-X(3))*X(4))

(k(1)*X(1)*(4.0626-X(4))-k(2)*X(3)*(4.0626-X(4))+k(3)*X(2)*(4.0626-X(4))-k(4)*X(3)*X(4)+k(5)*X(3)*(4.0626-X(4))-k(6)*(0.638-X(1)-X(2)-X(3))*X(4))];

%--------------------------------------------------

function f = jlzengObjFmincon(k,X0,Xexp)

%建立使用fmincon()进行参数估计的函数

tspan = [0 0.5 1 1.5 2 3 4 5 7 10];

[t X] = ode23s(@jlzengkinetic,tspan,X0,[],k);

f = sum((X(:,1)-Xexp(:,1)).^2) + sum((X(:,2)-Xexp(:,2)).^2)

...

+ sum((X(:,3)-Xexp(:,3)).^2) + sum((X(:,4)-Xexp(:,4)).^2);

%给物质浓度的计算值

%------------------------------------------------------

function f = jlzengObjLsqnonlin(k,X0,Xexp)

%在前面已经用fmincon求得参数的基础上使用Lsqnonlin对参数进行更精确的计算

tspan =[0 0.5 1 1.5 2 3 4 5 7 10];

[t X] = ode23s(@jlzengkinetic,tspan,X0,[],k);

f1 = X(:,1) - Xexp(:,1);

f2 = X(:,2) - Xexp(:,2);

f3 = X(:,3) - Xexp(:,3);

f4 = X(:,4) - Xexp(:,4);

f = [f1; f2; f3; f4];

%-------------------------------------------------------

运行结果发现1)两种方法求出的k值差别极大；2)不管将那组k值带入[t X] = ode23s(@jlzengkinetic,tspan,X0,[],k)求出的浓度与实际浓度差别很大；3)任意给定的k值的起始值对结果影响极大。

请高手帮忙！

如果还有没表述清楚的，请提出来，也可参看附件中的问题说明。