| 查看: 699 | 回复: 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求解微分方程组，然后用非线性最小二乘模拟参数，急求回复，谢谢大家！ |
回复此楼» 猜你喜欢
反铁磁体中的磁性切换：两种不同的机制已成功可视化
已经有0人回复
-
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有217人回复
-
新西兰Robinson研究所招收全奖PhD
已经有0人回复
-
石墨烯转移--二氧化硅衬底石墨烯
已经有0人回复
-
笼目材料中量子自旋液体基态的证据
已经有0人回复
-
数学教学论硕士可以读数学物理博士吗？
已经有0人回复
-
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生
已经有4人回复
-
澳门大学 应用物理及材料工程研究院 潘晖教授课题组诚招博士后
已经有11人回复
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