24小时热门版块排行榜    

Znn3bq.jpeg
查看: 1228  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 310求调剂 +6 666真好 2026-04-11 7/350 2026-04-11 20:49 by may_新宇
[考研] 305求调剂 +6 玛卡巴卡boom 2026-04-11 6/300 2026-04-11 19:00 by wutongshun
[考研] 085501机械专硕 302分 不挑专业求调剂 +7 汪某. 2026-04-09 7/350 2026-04-11 14:37 by luhong1990
[基金申请] 山东省基金2026 +4 jerry681 2026-04-08 5/250 2026-04-11 13:59 by laobibibi
[考研] 304求调剂(085602,过四级,一志愿985) +27 化工人999 2026-04-04 28/1400 2026-04-11 13:02 by 小晏唯
[考研] 求调剂 +6 archer.. 2026-04-09 8/400 2026-04-11 10:55 by zhq0425
[考研] 农学0904 312求调剂 +6 Say Never 2026-04-10 6/300 2026-04-11 10:33 by wwj2530616
[考研] 296求调剂 +13 汪!?! 2026-04-10 15/750 2026-04-11 10:31 by 逆水乘风
[考研] 生物学调剂 +8 小冉要努力 2026-04-10 9/450 2026-04-11 10:22 by wwj2530616
[考研] 调剂 +12 卷卷卷心菜_ 2026-04-09 13/650 2026-04-10 22:36 by Ftglcn90
[考研] 本科西工大 324求调剂 +4 wysyjs25 2026-04-10 4/200 2026-04-10 20:00 by 来看流星雨10
[考研] 一志愿华东师范生物学326分,求调剂 +8 刘墨墨 2026-04-09 8/400 2026-04-10 12:00 by pengliang8036
[考研] 材料复试求调剂 +20 xhhdjdjsjks 2026-04-09 20/1000 2026-04-10 10:25 by 孙小小12457
[考研] 一志愿华中农微生物,288分,三年实验经历 +10 代fish 2026-04-09 10/500 2026-04-10 09:49 by potato妹
[考研] 08工学 309分求调剂 +6 Yin DY 2026-04-08 6/300 2026-04-10 09:18 by Delta2012
[考研] 材料299专硕求调剂 +10 +21 2026-04-09 10/500 2026-04-09 17:34 by 1753564080
[考研] 考研求调剂 +4 雯??? 2026-04-08 4/200 2026-04-08 21:44 by 土木硕士招生
[考研] 0703调剂,一志愿天津大学319分 +23 haaaabcd 2026-04-05 26/1300 2026-04-08 16:19 by luoyongfeng
[考研] 0854求调剂 +9 亨氏番茄沙司 2026-04-06 10/500 2026-04-07 14:37 by shdgaomin
[考研] 第一志愿东南大学物理313,有科研竞赛获奖经历,希望物理复试调剂 +3 马内橙 2026-04-05 3/150 2026-04-06 10:32 by 蓝云思雨
信息提示
请填处理意见