| 查看: 747 | 回复: 1 | |||
| 【悬赏金币】回答本帖问题，作者szy242424将赠送您 5 个金币 | |||
[求助]
matalb软件 ode23s，非线性最小二乘法模拟动力学参数 已有1人参与
|
|||
|
function kk1 k0=[236,18,45.6,10,1.2,0.001,0.1,0.001,0.01,0.001,0.1,0.1]; lb=[236,18,45.6,10, 1.2, 0.001, 0.1,0.001,0.01,0.001,0.1,0.1]; ub=[inf,inf,inf,inf, 10.4,1.2, 2.4,1.2,0.8,0.1,4.5,4.5]; data=... [0 0.562205 0 0 34.2775 6 0.618817941 0 0 33.9805 12 0.797936454 0 0 31.941 18 1.554008384 0.9935 0 28.739 24 2.141789344 1.3815 0 26.3835 30 2.543955264 1.5745 1.4795 23.8955 36 3.017273616 1.908 1.7625 21.334 42 3.295696176 2.9885 2.038 19.128 60 3.274041088 4.693 3.262 11.2065 66 3.4070652 5.1295 3.581 8.9005 72 3.753546608 5.6395 4.041 6.4395 84 3.595773824 6.6735 4.71 1.496 ]; x0=data(1,2:end); tspan =[0,6,12,18,24,30,36,42,60,66,72,84]; yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) ]; ts=data(1:end,1); [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('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7 = %.9f ± %.9f\n',k(7),ci(7,2)-k(7)) fprintf('\tk8= %.9f ± %.9f\n',k(8),ci(8,2)-k(8)) fprintf('\tk9 = %.9f ± %.9f\n',k(9),ci(9,2)-k(9)) fprintf('\tk10 = %.9f ± %.9f\n',k(10),ci(10,2)-k(10)) fprintf('\tk11= %.9f ± %.9f\n',k(11),ci(11,2)-k(11)) fprintf('\tk12 = %.9f ± %.9f\n',k(12),ci(12,2)-k(12)) fprintf(' The sum of the squares is: %.9e\n\n',resnorm) [ts, ys] = ode23s(@KineticsEqs,ts,x0,[],k); yy = [data(:,2) data(:,3) data(:,4) data(:,5) ]; plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo'); hold on plot(ts,ys(:,2),'r',tspan,yy(:,2),'r*'); plot(ts,ys(:,3),'k',tspan,yy(:,3),'k+') legend('DCW的计算值','DCW的实验值','SA的计算值','SA的实验值','AA的计算值','AA的实验值','X的计算值','X的实验值') function f = ObjFunc(k,tspan,x0,yexp) % 目标函数 [~, Xsim] = ode23s(@KineticsEqs,tspan,x0,[],k); Xsim1=Xsim(:,1);%提取Xsim的第一列 Xsim2=Xsim(:,2);%提取Xsim的第二列 Xsim3=Xsim(:,3);%提取Xsim的第三列 Xsim4=Xsim(:,4);%提取Xsim的第四列 ysim(:,1) = Xsim1(2:end);%微分方程的第一个变量，从第二个解到最后一个解赋值给ysim的第一列 ysim(:,2) = Xsim2(2:end);%微分方程的第二个变量，从第二个解到最后一个解赋值给ysim的第二列 ysim(:,3) = Xsim3(2:end);%微分方程的第三个变量，从第二个解到最后一个解赋值给ysim的第三列 ysim(:,4) = Xsim4(2:end);%微分方程的第四个变量，从第二个解到最后一个解赋值给ysim的第四列 f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,4)-yexp(:,4))]; function dCdt = KineticsEqs(~,C,k) % ODE模型方程，C1、C2、C3、C4分别为底物浓度、产物丁二酸浓度、产物乙酸浓度、DCW值，t为导数 dC1dt =(0.1933862*C(4)/(C(4)+0.78+C(4)^2/296.9)*(1-C(2)/k(1))^k(2)*(1-C(3)/k(3))^k(4))*C(1); dC2dt =k(5)*dC1dt+k(6)*C(1);%第二个微分方程，等号前面是C(2)对t的微分 dC3dt =k(7)*dC1dt+k(8)*C(1);%第三个微分方程，等号前面是C(3)对t的微分 dC4dt =-((1/k(9))*dC1dt+k(10)*C(1)+(1/k(11))*dC2dt+(1/k(12))*dC3dt); dCdt = [dC1dt; dC2dt;dC3dt;dC4dt];%微分方程组. 各位大神，帮我看一下这个代码有什么错误没有，目标是用ode23s求解微分方程组，然后用非线性最小二乘模拟参数，急求回复，谢谢大家！ |
回复此楼» 猜你喜欢
散金币，求好运，祝面上顺利！
已经有80人回复
-
华东师范大学芯片设计徐珑真实水平咋样
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有190人回复
-
调剂
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有0人回复
-
中德博士后2026项目5月30日截止- 国内博士后入站、德国两年全职
已经有0人回复
-
面向全球申请人的中德博士后2026项目5月30日截止- 国内博后入站德国两年全职
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有2人回复
-
重庆交大材料学院26年硕士招生复试实施细则已出。欢迎加入光子学微结构与器件课题组。
已经有0人回复
-
北工大宋晓艳教授团队招收2026年9月入学博士【学博和专博名额均有】！
已经有1人回复
-
天津城建大学理学院化学、环境工程、材料物理等方向接收10人
已经有0人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
感谢参与，应助指数 +1
|
参考下面1stOpt计算的结果： Root of Mean Square Error (RMSE): 0.673324275669162 Sum of Squared Residual: 19.9480855290377 Correlation Coef. (R): 0.983872215020488 R-Square: 0.968004535489321 Parameter Best Estimate -------------------- ------------- k1 2255.98568602798 k2 2006.16856712862 k3 123314855.293705 k4 11.7085151403374 k5 1.20000000001939 k6 0.00942914835572298 k7 0.100000000023396 k8 0.0201718372859985 k9 0.799999999987507 k10 0.0999392603337909 k11 4.49999420593383 k12 0.75747434507548 |
2楼2021-03-10 10:18:00
