A+B=C  A+W=R  R+B=C+W
先考虑的三个反应都是可逆反应，有以下几个式子：编写程序如下
%30℃
clear all
clc
global keq1 keq2 t0 c0

keq1 = 3.94;         % 平衡常数k1
keq3 = 1.08;   % 平衡常数k3

t0 = [5,10, 15, 30, 45, 60,90,120,160,200,240,300,360];
c0 = [5.939601,5.203973,0,0.022484,0.078956]        % A B C R W 组成的初始浓度
ci=[3.742508811        4.213566482        0.990406938        0.252498405        0.090415512;
3.045716838        3.387013322        1.816960097        0.215546741        0.081072022;
2.504754523        2.901991595        2.301981824        0.193427102        0.099462604;
1.65098702        2.036686653        3.167286766        0.170236582        0.095144044;
1.271633833        1.655460403        3.548513016        0.124730781        0.10334349;
1.129887499        1.630815394        3.573158025        0.129798695        0.097764543;
1.01239208        1.460748979        3.74322444        0.144649391        0.073614958;
1.007417603        1.429180292        3.774793127        0.123477124        0.094249307;
0.990005723        1.499216939        3.704756481        0.135524067        0.076072022;
0.970073726        1.493966669        3.71000675        0.130791851        0.076387812;
0.96903024        1.429705429        3.77426799        0.119510017        0.101202216;
0.980873385        1.457988124        3.745985296        0.124189274        0.097022161;
0.98477262        1.453402965        3.750570455        0.1200936         0.0762608];

% cA cC cR 对应时间数值

k0=[1 1 1];   % k1+, k2, k3+ 初始向值
lb = [0 0 0];
ub = [1000  1000 500];  % 上下限

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag] = lsqnonlin(@ObjFunc,k0,lb,ub,[],ci);      

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

% ------------------------------------------------------------------
function f = ObjFunc(k,ci)    % 目标函数
global t0  c0

[t,c_cal] = ode45(@Euqations,t0,c0,[],k);
f= sum((c_cal-ci)^2);

% ------------------------------------------------------------------
function dcdt = Euqations(t,c,k)
global keq1 keq2

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

cA=c(:,1);
cB=c(:,2);
cC=c(:,3);
cR=c(:,4);
cW=c(:,5);

dcAdt = k1plus*cA*cB-k1minus*cC+k2*cA*cW;
dcCdt = -k1plus*cA*cB+k1minus*cC+k3plus*cB*cR-k3minus*cC*cW;
dcRdt = k2*cA*cW-k3plus*cB*cR+k3minus*cC*cW;
dcdt = [dcAdt; dcCdt;dcRdt];
但是k的初值一直影响结果，程序要怎么改一下初值对结果影响不大（这个程序有问题，） 返回小木虫查看更多