24小时热门版块排行榜    

查看: 1200  |  回复: 4

anren

铜虫 (初入文坛)

[求助] 求助matlab拟合动力学K值 已有1人参与

动力学模型为
dC(1)/dt=    -((k(1)+k(2)))*C(1)
dC(2)/dt=   k(1)*C(1)-k(3)*C(2)-k(4)*C(2)
dC(3)/dt=    k(2)*C(1)-k(5)*C(3)
dC(4)/dt=    -k(6)*C(4)+k(3)*C(2)
dC(5)/dt=     k(4)*C(2)+k(5)*C(3)-k(7)*C(5)
dC(6)/dt=     k(6)*C(4)+k(7)*C(5) ;

实验数据   
%  t      C(1)      C(2)    C(3)   C(4)    C(5)     C(6) %   
   ExpData =[
         0       238.5300         0         0           0         0            0                                               
   22.0000   77.1500    8.3600   13.9900        0    0.5100         0
   36.0000   58.4500   13.0200   27.5800    0.0300    0.9200         0
   55.0000   33.0900   16.8200   48.3100    0.0500    1.7100    0.0100
   76.0000   15.9600   25.5900   55.2300    0.1200    3.0500    0.0500
   96.0000    1.5500   25.7900   63.3900    0.4200    8.8400         0
  115.0000    0.0400   19.1900   46.9500    1.2800   32.4400    0.1000
  134.0000         0   17.8700   25.0700    2.0000   54.8500    0.2100
  154.0000    0.0200    6.2800   15.7100    2.5300   74.9700    0.4900
  173.0000         0    1.6900    3.0300    2.9600   91.3100    0.9900
  177.0000         0    1.1000    1.3300    2.9300   93.6800    0.9500
  178.0000         0    0.5400    0.3600    3.1000   94.9000    1.0800
  188.0000         0    0.3700    0.5000    3.0400   94.7400    1.3500 ]

哪位大神能帮忙拟合一下k值呢?为什么我一直得不出结果
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

youmengtian

铜虫 (初入文坛)

【答案】应助回帖

感谢参与,应助指数 +1
你这个是怎么思路做的?我这里有一个简单的程序,你可以看一眼。
function k1k2k3
format long
clear all
clc
tspan = [0  30  50  80 140 170 200 230 260 290 320 360 400 460 520 580 600];
x0 = [9; 0];
k0 = [0.001  0  0.08];   
lb = [0  0  0];
ub = [1  1  1];

data=...
    [

30        7.939        1.458
50        7.687        1.535
80        7.289        1.602
140        6.658        1.717
170        6.531        1.722
200        6.218        1.671
230        5.979        1.620
260        5.591        1.550
290        5.414        1.488
320        4.968        1.433
360        4.692        1.350
400        4.438        1.319
460        4.144        1.294
520        4.041        1.294
580        4.052        1.287
600        4.030        1.292
];
yexp = data(:,2:3);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)

function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);

%size(ysim(:,1));
%size(ysim(:,2));
%size(yexp(:,1));
%size(yexp(:,2));

f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2))];

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
dCAdt = -k(1)*C(1)+k(2)*C(2);
dCBdt = k(1)*C(1)-k(2)*C(2)-k(3)*C(2);
dCdt = [dCAdt; dCBdt];
2楼2015-06-26 19:41:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

matlab编程

禁虫 (小有名气)

本帖内容被屏蔽

3楼2015-06-27 12:34:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anren

铜虫 (初入文坛)

引用回帖:
2楼: Originally posted by youmengtian at 2015-06-26 19:41:28
你这个是怎么思路做的?我这里有一个简单的程序,你可以看一眼。
function k1k2k3
format long
clear all
clc
tspan = ;
x0 = ;
k0 = ;   
lb = ;
ub = ;

data=...
    ;
yexp = data(:,2:3);

=  ...

这个是我照书上程序改的,可是计算不出结果。
function KineticsEst5
clear all
clc
k0 = [0.5  0.5  0.5  0.5  0.5 0.5 0.5];         % 参数初值
lb = [0  0  0  0  0 0 0];                   % 参数下限
ub = [+inf  +inf  +inf  +inf  +inf +inf  +inf];    % 参数上限
C0 = [ 238.53    0   0   0   0   0      ];
ExpData =[0 238.53  0   0   0   0   0
22  77.15   8.36    13.99   0   0.51    0
36  58.45   13.02   27.58   0.03    0.92    0
55  33.09   16.82   48.31   0.05    1.71    0.01
76  15.96   25.59   55.23   0.12    3.05    0.05
96  1.55    25.79   63.39   0.42    8.84    0
115 0.04    19.19   46.95   1.28    32.44   0.1
134 0   17.87   25.07   2   54.85   0.21 ]
t= ExpData(:,1)'
yexp = ExpData(:,2:7);                  % yexp: 实验数据[x1        x4        x5        x6]

% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],C0,yexp);
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('\tk5 = %.4f\n',k(6))
fprintf('\tk5 = %.4f\n',k(7))

fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],C0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
Output

% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],C0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
Output


% ------------------------------------------------------------------
function f = ObjFunc4Fmincon(k,C0,yexp)
tspan = [0    22    36    55    76    96   115   134 ];
[t C] = ode45(@KineticEqs,tspan,C0,[],k);   
y(:,1) = C(:,1);
y(:,2:6) = C(:,2:6);
f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)   ...
    + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2) ...
    + sum((y(:,5)-yexp(:,5)).^2) + sum((y(:,6)-yexp(:,6)).^2);

% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,C0,yexp)
tspan = [0    22    36    55    76    96   115   134 ];
[t C] = ode45(@KineticEqs,tspan,C0,[],k);   
y(:,1) = C(:,1);
y(:,2:6) = C(:,2:6);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f6 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4;f5;f6];

% ------------------------------------------------------------------
function dCdt = KineticEqs(t,C,k)
dCdt =  ...
       [  ( -((k(1)+k(2)))*C(1) )
           (k(1)*C(1)-k(3)*C(2)-k(4)*C(2) )
           ( k(2)*C(1)-k(5)*C(3) )
          ( -k(6)*C(4)+k(3)*C(2) )
          ( k(4)*C(2)+k(5)*C(3)-k(7)*C(5) )
        ( k(6)*C(4)+k(7)*C(5) )];
4楼2015-06-29 10:07:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anren

铜虫 (初入文坛)

引用回帖:
3楼: Originally posted by matlab编程 at 2015-06-27 12:34:16
这个我才做过。

那快点把程序分享一下啊
5楼2015-06-29 10:10:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 anren 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 307求调剂 +3 73372112 2026-02-28 4/200 2026-02-28 23:49 by 迷糊CCPs
[考研] 272求调剂 +3 材紫有化 2026-02-28 3/150 2026-02-28 22:52 by ms629
[考研] 304求调剂 +3 52hz~~ 2026-02-28 4/200 2026-02-28 21:41 by gaoxiaoniuma
[考研] 290求调剂 +5 材料专硕调剂; 2026-02-28 6/300 2026-02-28 21:40 by gaoxiaoniuma
[考博] 26申博 +4 想申博! 2026-02-26 4/200 2026-02-28 21:37 by limorning
[考研] 264求调剂 +3 巴拉巴拉根556 2026-02-28 3/150 2026-02-28 21:31 by gaoxiaoniuma
[考研] 311求调剂 +8 南迦720 2026-02-28 8/400 2026-02-28 21:30 by gaoxiaoniuma
[考研] 284求调剂 +4 天下熯 2026-02-28 4/200 2026-02-28 21:13 by gaoxiaoniuma
[考研] 298求调剂 +8 人间唯你是清欢 2026-02-28 11/550 2026-02-28 20:26 by L135790
[考研] 276求调剂 +3 路lyh123 2026-02-28 4/200 2026-02-28 19:45 by 路lyh123
[教师之家] 版面费该交吗 +15 苹果在哪里 2026-02-22 18/900 2026-02-28 18:20 by mibaomingg
[考研] 材料调剂 +3 爱擦汗的可乐冰 2026-02-28 3/150 2026-02-28 18:06 by houyaoxu
[高分子] 求环氧树脂研发1名 +3 孙xc 2026-02-25 11/550 2026-02-28 16:57 by ichall
[考研] 265分求调剂不调专业和学校有行学上就 +4 礼堂丁真258 2026-02-28 6/300 2026-02-28 16:18 by 求调剂zz
[考研] 寻找调剂 +3 LYidhsjabdj 2026-02-28 3/150 2026-02-28 12:59 by miniwendy
[硕博家园] 博士自荐 +6 科研狗111 2026-02-26 9/450 2026-02-28 12:32 by seaskyy
[考研] 272求调剂 +3 田智友 2026-02-28 3/150 2026-02-28 12:31 by 王加浩to
[考研] 298求调剂 +4 axyz3 2026-02-28 4/200 2026-02-28 11:21 by wang_dand
[基金申请] 面上可以超过30页吧? +12 阿拉贡aragon 2026-02-22 13/650 2026-02-26 22:09 by Hahaxia
[硕博家园] 【博士招生】太原理工大学2026化工博士 +4 N1ce_try 2026-02-24 8/400 2026-02-26 08:40 by N1ce_try
信息提示
请填处理意见