| 查看: 910 | 回复: 3 | ||
| 【悬赏金币】回答本帖问题，作者自制奶酪将赠送您 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;]; |
回复此楼» 猜你喜欢
中国科学院东莞材料科学与技术研究所-2026年博士招生-吴昊研究员-磁学与自旋电子学
已经有0人回复
-
《电磁学》教材推荐
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有100人回复
-
【急招】合肥工大核聚变材料计算方向2026级工程博士生
已经有4人回复
-
大豆异黄酮分离
已经有0人回复
-
湖南大学材料学院急招2026年博士生，临时增加一名博士联培指标
已经有10人回复
-
天津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生1-2名
已经有1人回复
-
中国科学院物理研究所谌志国研究员团队招收2027年博士研究生
已经有3人回复
-
2026年中德博士后交流项目 - 新型量子和磁性材料：材料制备表征和中子散射研究
已经有12人回复
hzlhm
至尊木虫 (著名写手)
- 应助: 387 (硕士)
- 金币: 17823.5
- 红花: 53
- 帖子: 2879
- 在线: 606.3小时
- 虫号: 1517335
- 注册: 2011-11-30
- 性别: GG
- 专业: 常微分方程与动力系统
2楼2021-01-11 19:19:58
3楼2021-04-24 21:45:22
4楼2021-04-24 21:48:28
10