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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 理论与计算方法 » 用matlab求动力学参数k1-k7,金币没有，太寒碜请见谅
41/1返回列表
查看: 831  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
【悬赏金币】回答本帖问题，作者自制奶酪将赠送您 30 个金币

自制奶酪

铜虫 (小有名气)

[求助] 用matlab求动力学参数k1-k7,金币没有，太寒碜请见谅 已有1人参与

希望木虫大神帮忙用Matlab解解方程，可以酬金支付，

t/h      C1        C2       C3       C4        C5
0       100       0      0      0      0
0.5      18.7     48.9     19.3      9.7       0.2
1       6.5       31.3     34.7     19.6      2.6
1.5      2.7       15.3     41.2     28.7      5.7
2       2.1       6.9      40.9      34.5      9.3
3       1.7       2.2      37.2      38.4     17.1
4       1.5       1.8      31.7      38.2      24.8
5       1.2       1.7      26.1      36.3      30.5
6       0.2       0.5      19.4      33.4      36.7

微分方程组

dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C6;
dC4dt = k3*C3-k4*C4;
dC5dt = k4*C4;
dC6dt = k5*C2+k6*C3+k7*C1;
我网上依葫芦画瓢填的报错结果出不来

代码：

function odes_fit
format long
clear all
clc


k0 = [0 0 0 0 0 0];   
lb = -[1 1 1 1 1 1]*1e9;
ub = [1 1 1 1 1 1]*1e9;

data=...
    [0         100             0            0            0            0            0
         0.5    18.7   48.9   19.3    9.7   0.2   0.0;
          1    6.5   31.3   34.7   19.6   2.6   0.0;
          1.5    2.7    15.3   41.2   28.7   5.7   0.05;
          2    2.1    6.9    40.9   34.5   9.3   0.1;
          3    1.7    2.2    37.2   38.4   17.1   0.12;
          4    1.5    1.8    31.7   38.2   24.8   0.16;
          5    1.2    1.7    26.1   36.3   30.5   0.22;
          6    0.2    0.5    19.4   33.4   36.7   0.3;
];
x0=data(1,2:end);
tspan = [data(:,1)'];
yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) data(2:end,6) data(2:end,7)];

[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 \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk5 = %.9f \n',k(6))
fprintf('\tk5 = %.9f \n',k(7))


figure(1)

ts=0(max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4) data(:,5) data(:,6) data(:,7)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo');
figure(2)
plot(ts,ys(:,2),'r',tspan,yy(:,2),'ro');
figure(3)
plot(ts,ys(:,3),'k',tspan,yy(:,3),'ko');
figure(4)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'go');
figure(5)
plot(ts,ys(:,5),'m',tspan,yy(:,5),'mo');
figure(6)
plot(ts,ys(:,6),'h',tspan,yy(:,6),'ho');
figure(7)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro',ts,ys(:,3),'k',tspan,yy(:,3),'ko',ts,ys(:,4),'g',tspan,yy(:,4),'go',ts,ys(:,5),'m',tspan,yy(:,5),'mo',ts,ys(:,5),'h',tspan,yy(:,5),'ho'),
legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值','C5的计算值','C5的实验值','C6的计算值','C6的实验值','Location','best');



function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
Xsim5=Xsim(:,5);
Xsim6=Xsim(:,6);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
ysim(:,6) = Xsim6(2:end);



f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,5)-yexp(:,5))  (ysim(:,6)-yexp(:,6)) ...
    (ysim(:,6)-yexp(:,6))];

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
C1=C(1);C2=C(2);C3=C(3);C4=C(4);
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);
dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C6;
dC4dt = k3*C3-k4*C4;
dC5dt = k4*C4;
dC6dt = k5*C2+k6*C3+k7*C1;

dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt;];
回复此楼

» 猜你喜欢
简单快乐
1楼 2021-01-06 13:34:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
独孤神宇: 金币+10, 鼓励交流 2021-01-11 20:05:59
如方程没有错误的话，其k值为
k1=2.9355； k2=1.318； k3=0.44525； k4=0.2476； k5=2.2315e-14； k6=2.2214e-14； k7=0.47562
C1的拟合精度R^2=0.99696
C2的拟合精度R^2=0.98746
C3的拟合精度R^2=0.74181
C4的拟合精度R^2=0.71926
C5的拟合精度R^2=0.94598
赞一下(2人) 回复此楼
QQ：2120156492
2楼2021-01-11 19:19:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

梅子炒米粉

新虫 (小有名气)
请问楼主问题解决了吗？方便交流交流吗？

发自小木虫Android客户端
赞一下 回复此楼
3楼2021-04-24 21:45:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

梅子炒米粉

新虫 (小有名气)
引用回帖:
2楼: Originally posted by hzlhm at 2021-01-11 19:19:58
如方程没有错误的话，其k值为
k1=2.9355； k2=1.318； k3=0.44525； k4=0.2476； k5=2.2315e-14； k6=2.2214e-14； k7=0.47562
C1的拟合精度R^2=0.99696
C2的拟合精度R^2=0.98746
C3的拟合精度R^2=0.74181
C4的 ...

请问前辈知道用matlab该怎么改上面的模型嘛？

发自小木虫Android客户端
赞一下 回复此楼
4楼2021-04-24 21:48:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 自制奶酪 的主题更新
41/1返回列表
不应助 确定回帖应助 (注意：应助才可能被奖励，但不允许灌水，必须填写15个字符以上)
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭