matlab拟合不稳定初值影响大,matlab拟合方程参数时初值的选择 - 计算模拟 - 小木虫 - 学术 科研 互动社区...

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的初值一直影响结果，程序要怎么改一下初值对结果影响不大(这个程序有问题，)