应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】反应动力学求解
141/2返回列表12下一页
查看: 1820  |  回复: 13
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

jlzeng

木虫 (正式写手)

[交流] 【求助】反应动力学求解

我在求反应动力学时遇到一个非常棘手的问题，请大家帮忙解决一下，先对每一个关注本贴的人表示感谢！

问题描述：
       有三个可逆反应组成的反应体系如下：
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

我用MATLAB求解某一条件下的动力学常数k，程序如下：

clear all
k0 = [1 1 3 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值的起始值对结果影响极大。
请高手帮忙！
如果还有没表述清楚的，请提出来.

[ Last edited by sunxiao on 2009-3-8 at 11:17 ]
回复此楼

» 收录本帖的淘帖专辑推荐

matlab

» 猜你喜欢
沉默坚守
1楼 2009-03-04 10:29:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jlzeng

木虫 (正式写手)
★
kuhailangyu(金币+1,VIP+0):帮你加亮显示了，话题不错，讲的也挺详细的，耐心等待相关高手来解答吧！ 3-4 11:02
猜想是刚性方程在由OED求数值解时出现了问题（多凸还是不连续时只在局部求解而不是对全局求最优解），不知道能不能强制在0-50之间运行每一个可能的解。
我是学生物的，现在搞化工，对计算实在是一窍不通啊！这一部分的内容是论文里必须的，还望各位高手出手相助，非常感谢！！
赞一下(10人) 回复此楼
沉默坚守
2楼2009-03-04 10:45:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jlzeng

木虫 (正式写手)
怎么没有人回复啊？
~~~~(>_<~~~~
赞一下 回复此楼
沉默坚守
3楼2009-03-05 09:10:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hitzhang

木虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+3,VIP+0):欢迎交流，得到楼主认可后，将继续追加奖励 3-5 23:33
jlzeng(金币+10,VIP+0):很有用，但还没有完全解决，非常感谢！ 3-6 15:48
我看了一下你的方法，感觉很新颖，是一种纯粹的数值解法，想必你的MATLAB功底不浅。你主要是想通过两个优化工具箱函数利用最小二乘法解微分方程组进而求出速率常数。但是效果不好，我想你下面的分析是有道理的，就算法有效性而言fmincon这个函数口碑不是很好，多数都是用LINGO软件求优化问题，更别说求微分方程组了，而且透明度也不高好像在解决一个黑箱问题，求解对初始条件很敏感。

我这里有个想法仅供你参考：

这个问题的最终目的是求最佳速率常数以尽量满足微分方程组，现在，你有了各个物种随时间变化的浓度，如果知道浓度随时间的变化速率，就可以直接代入微分方程组应用最小而乘求出速率常数k。
浓度变化速率的最佳估计：可以根据你的实验数据构造出在实验范围内符合较好的各个物种的浓度方程。以下是拟和结果：

///////////////////////////////////////////////////////////////////////////////
TG:
       f(x) = a*exp(b*x) + c*exp(d*x)
Coefficients (with 95% confidence bounds):
       a =      0.5502  (0.5149, 0.5856)
       b =      -3.983  (-4.797, -3.17)
       c =     0.08751  (0.05922, 0.1158)
       d =     -0.1857  (-0.2901, -0.08124)

Goodness of fit:
  SSE: 0.0004627
  R-square: 0.9985
  Adjusted R-square: 0.9978
  RMSE: 0.008782



DG:
       f(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2)
Coefficients (with 95% confidence bounds):
       p1 =     0.01908  (0.01341, 0.02475)
       p2 =     0.05717  (0.01411, 0.1002)
       p3 =  -1.43e-007  (-0.0008759, 0.0008756)
       q1 =     -0.2444  (-0.8719, 0.3831)
       q2 =      0.1788  (0.04132, 0.3162)

Goodness of fit:
  SSE: 1.816e-005
  R-square: 0.9979
  Adjusted R-square: 0.9963
  RMSE: 0.001906


MG:
       f(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2)
Coefficients (with 95% confidence bounds):
       p1 =     0.07364  (0.04582, 0.1015)
       p2 =      0.3095  (-0.1137, 0.7328)
       p3 =  1.802e-005  (-0.00794, 0.007976)
       q1 =      0.9642  (-1.498, 3.426)
       q2 =      0.6633  (0.1897, 1.137)

Goodness of fit:
  SSE: 0.0001089
  R-square: 0.9933
  Adjusted R-square: 0.988
  RMSE: 0.004666


ME:
       f(x) = (p1*x + p2) / (x^2 + q1*x + q2)
Coefficients (with 95% confidence bounds):
       p1 =        2169  (-1.654e+004, 2.088e+004)
       p2 =      0.9198  (-23.22, 25.06)
       q1 =        1265  (-9708, 1.224e+004)
       q2 =       406.4  (-3050, 3862)

Goodness of fit:
  SSE: 0.003249
  R-square: 0.9986
  Adjusted R-square: 0.9978
  RMSE: 0.02327

///////////////////////////////////////////////////////////////////////

这几个浓度方程的拟和结果如如附件中所示。

然后就可以计算在实验时间点上TG, DG, MG, ME的浓度对时间的导数了，结果如下：////////////////////////////////////////////////////////////////////////

>> dTGdt=[-2.20812330989885;-0.313903194394831;-0.0543072642775238;-0.0178664548644437;-0.0119669262873315;-0.00932216291313879;-0.00773110042346866;-0.00642090783831632;-0.00442930197092527;-0.00253772213004114]'

dTGdt =

   -2.2081   -0.3139   -0.0543   -0.0179   -0.0120   -0.0093   -0.0077   -0.0064   -0.0044   -0.0025

>> dDGdt=[0.319813865326457;-0.0194568056055164;-0.0513088546928754;-0.0279068843505723;-0.0164124937537988;-0.00737171342598024;-0.00411902140843388;-0.00261612659035434;-0.00131890646666286;-0.000639095343223927;]'

dDGdt =

    0.3198   -0.0195   -0.0513   -0.0279   -0.0164   -0.0074   -0.0041   -0.0026   -0.0013   -0.0006

>> dMGdt=[0.466367050839715;0.0999072143408668;0.00933499024234685;-0.00972884650823934;-0.0127379996363272;-0.0104546020606972;-0.00764618068387975;-0.00566858666344666;-0.00339187374122912;-0.00186232961806277;]'

dMGdt =

    0.4664    0.0999    0.0093   -0.0097   -0.0127   -0.0105   -0.0076   -0.0057   -0.0034   -0.0019

>> dMEdt=[5.28500993498458;0.818198974835815;0.315950868952596;0.165592835535058;0.101261692838514;0.0484953510802743;0.0278523652553689;0.0177021334686043;0.00842980322026890;0.00327667396220621;]'

dMEdt =

    5.2850    0.8182    0.3160    0.1656    0.1013    0.0485    0.0279    0.0177    0.0084    0.0033

>>

////////////////////////////////////////////////////////////////////////////

之后的工作就是把这四组数据和实验数据代入微分方程组，求出最佳速率常数。
赞一下(20人) 回复此楼
4楼2009-03-05 21:16:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hitzhang

木虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+2,VIP+0):欢迎交流，得到楼主认可后，将继续追加奖励 3-5 23:33
jlzeng(金币+15,VIP+0):这个问题得到hitzhang很多启示,非常感谢！！ 4-22 10:23
下面是对速率常数k的求解：

你那个微分方程组可以写成下面这种矩阵形式：
C'=A*k

下面的人物就是如何解出k
按照楼上所说的方法：
C'=

   -2.2081
    0.3198
    0.4664
    5.2850
   -0.3139
   -0.0195
    0.0999
    0.8182
   -0.0543
   -0.0513
    0.0093
    0.3160
   -0.0179
   -0.0279
   -0.0097
    0.1656
   -0.0120
   -0.0164
   -0.0127
    0.1013
   -0.0093
   -0.0074
   -0.0105
    0.0485
   -0.0077
   -0.0041
   -0.0076
    0.0279
   -0.0064
   -0.0026
   -0.0057
    0.0177
   -0.0044
   -0.0013
   -0.0034
    0.0084
   -0.0025
   -0.0006
   -0.0019
    0.0033

A=


   -2.5919         0         0         0         0         0
    2.5919         0         0         0         0         0
         0         0         0         0         0         0
    2.5919         0         0         0         0         0
   -0.4516    0.0135         0         0         0         0
    0.4516   -0.0135   -0.3250    0.1338         0         0
         0         0    0.3250   -0.1338   -0.3716    0.0315
    0.4516   -0.0135    0.3250    0.1338    0.3716   -0.0315
   -0.2663    0.0119         0         0         0         0
    0.2663   -0.0119   -0.2275    0.1852         0         0
         0         0    0.2275   -0.1852   -0.4055    0.0459
    0.2663   -0.0119    0.2275    0.1852    0.4055   -0.0459
   -0.1925    0.0090         0         0         0         0
    0.1925   -0.0090   -0.1692    0.1983         0         0
         0         0    0.1692   -0.1983   -0.3831    0.0515
    0.1925   -0.0090    0.1692    0.1983    0.3831   -0.0515
   -0.1289    0.0069         0         0         0         0
    0.1289   -0.0069   -0.1271    0.2077         0         0
         0         0    0.1271   -0.2077   -0.3610    0.0558
    0.1289   -0.0069    0.1271    0.2077    0.3610   -0.0558
   -0.1145    0.0056         0         0         0         0
    0.1145   -0.0056   -0.1107    0.1956         0         0
         0         0    0.1107   -0.1956   -0.3256    0.0540
    0.1145   -0.0056    0.1107    0.1956    0.3256   -0.0540
   -0.0898    0.0043         0         0         0         0
    0.0898   -0.0043   -0.0873    0.1932         0         0
         0         0    0.0873   -0.1932   -0.3048    0.0545
    0.0898   -0.0043    0.0873    0.1932    0.3048   -0.0545
   -0.0847    0.0031         0         0         0         0
    0.0847   -0.0031   -0.0737    0.1660         0         0
         0         0    0.0737   -0.1660   -0.2516    0.0484
    0.0847   -0.0031    0.0737    0.1660    0.2516   -0.0484
   -0.0631    0.0029         0         0         0         0
    0.0631   -0.0029   -0.0681    0.1672         0         0
         0         0    0.0681   -0.1672   -0.2455    0.0490
    0.0631   -0.0029    0.0681    0.1672    0.2455   -0.0490
   -0.0569    0.0025         0         0         0         0
    0.0569   -0.0025   -0.0617    0.1625         0         0
         0         0    0.0617   -0.1625   -0.2344    0.0480
    0.0569   -0.0025    0.0617    0.1625    0.2344   -0.0480


k=A\C'


    1.0065
   17.2251
    0.6860
    0.1923
    0.9624
    5.9171
赞一下(16人) 回复此楼
5楼2009-03-05 22:01:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jlzeng

木虫 (正式写手)
非常感谢hitzhang的回复！让您费了不少心思~
您的解题方法我曾经也想试试，不过由于存在下面的问题，没有解决：
1）        我用MATLAB直接拟合的曲线不是太好，表现在根据拟合曲线求出的倒数不是趋近于0（应该趋近于0才对的）；hitzhang根据图形趋势选用合适的曲线类型进行拟合所求出的倒数非常好；
2）        我对此处用C'=A*k来求解线性方程的原理想不通。因为k在此处是由6个未知数组成的向量，而由C'=A*k组成的方程组有4*10=40个方程，用40个方程如何才能求出6个未知数而不出现无解的呢？

我把您所求出的k值带入MATLAB，用如下程序（见附录）求此条件下各物质浓度分布时出现“Warning: Failure at t=5.592464e-001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.986842e-015) at time t.
> In ode23s at 444”的警告。出现上述结果可能有两个原因：1）所求k值不对；2）所求k值没有问题，但是由于方程本身的特殊性，在调用ODE23s求解时，无法计算。到底是那一种原因，我暂时还没法确证。
在此我想先请hitzhang跟我讲讲用C'=A*k来求解线性方程的原理，如果原理说得通，我就用相同的方法把其他几组数据算出来，相互之间一比对就可确定该方法的正确性了！如果方便的话，我想请hitzhang把您的电话或手机号码发到我的邮箱 jlzeng@home.ipe.ac.cn ,在电话里交流可能更方便快捷些，非常感谢！

        附录：
X0=[0.638;0;0;0];                        %t=0时刻四种物质的初始浓度
  k=[1.0065   17.2251    0.6860    0.1923    0.9624    5.9171]
  tspan = [0 0.5 1 1.5 2 3 4 5 7 10]; %3 4 5 7 10 15 20 30 60
  [t X] = ode23s(@jlzengkinetic,tspan,X0,[],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))];

[ Last edited by jlzeng on 2009-3-6 at 15:47 ]
赞一下 回复此楼
沉默坚守
6楼2009-03-06 15:45:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hitzhang

木虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+8,VIP+0):重奖，欢迎来做我们的专家顾问 3-6 23:58
kuhailangyu(金币+0,VIP+0):你参与的主题已有新的进展，敬请关注，谢谢！ 3-14 15:57
微分方程组：
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]
有4个方程，可简写为：c'=a*k。你测定了10个时刻4个物种的浓度，可计算出每个时刻的矩阵a，另外这10个时刻各个物种的反应速率c'已估计出。在每个时刻上满足c'=a*k，有4个方程，那么十个时刻共有40个方程。可把这40个方程简写为：C'=A*K。数学上称这类方程为超定方程，即方程组中方程的个数多于未知量的个数，一般情况超定方程无解，只能找到一个最近似的解，主要应用最小二乘法，使方程左右两边尽量相等。具体算法可参阅相关文献。

在matlab中解这个方程显得异常简单，K=A\C'。

一个严重的问题：结果代入微分方程模型进行述职求解时和实验偏差较大。我想一个很重要的原因是我把矩阵A算错了，因为繁琐避免不了马虎，如果矩阵A没有错误，我想最可能的原因是微分方程模型不够合理。仅供参考。
赞一下(20人) 回复此楼
7楼2009-03-06 19:36:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jlzeng

木虫 (正式写手)
收到！
下周有个面试，还要准备ppt，等ppt和面试搞定了再回来做这个题，验证您的方法！
你的想法（微分方程模型可能不合理）是我以前没想到，但很有可能的！
非常感谢您的意见！
等我算完后再将结果告诉您
赞一下 回复此楼
沉默坚守
8楼2009-03-06 22:50:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jlzeng

木虫 (正式写手)
★
sunxiao(金币+1,VIP+0):鼓励一下，努力 3-15 02:16
这几天按照 hitzhang的方法，算了其他的几组数据，好像还是不行。
考虑到模型可能有问题，添加了一个可能的方程，结果还是有问题。
从图上可以看出，反应在开始阶段比较快，后面就比较慢了（这种现象在催化剂浓度更高、温度更高时，更加明显），这是不是意味着前面的方程和后面的方程不一样（例如前面是二级反应，后面是0级或一级反应）？这样就不能用一个方程表述了？也就是说我所建立的模型只适用于快速反应的二级方程阶段（即反应起始阶段）？后续阶段的数据不能用于模型求解？
赞一下(10人) 回复此楼
沉默坚守
9楼2009-03-14 15:54:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snipher950

木虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+3,VIP+0):谢谢参与，欢迎常来 3-15 08:48
jlzeng(金币+10,VIP+0):谢谢您的建议！ 3-16 10:57
jlzeng(金币+10,VIP+0):您让我对数值解法又增加了认识，非常感谢！ 4-22 10:24
粗略看了楼上两位的讨论，其实微分方程组的参数估计基本就是你们讨论的两种方法：积分法和微分法。
所谓积分法（楼主所用），就是根据预设的参数初值用积分的方法求解微分方程，然后与实验所得浓度值最小二乘拟合。这种方法与初值很有关系，当然同时更决定于你用哪种优化算法去拟合。目前给出的大多算法都是局部优化算法，自然会出现楼主所述情况。
所谓微分法（hitzhang所用），就是将所测的试验浓度曲线微分，然后把相关浓度数值和预设参数直接带入微分方程组，得到Dx/dt，再与浓度曲线微分结果拟合以求得参数，这种情况在实验数据充足的情况下，更加容易拟合，计算速度也更快一些。

参数估计遇到这样的过程都很正常，可以做灵敏度分析来选择灵敏参数进行估计，而将不灵敏的设为某个定值，便量少，便于分析一些。

仅供参考。
赞一下(21人) 回复此楼
10楼2009-03-15 08:33:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jlzeng 的主题更新
141/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭