| 查看: 760 | 回复: 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求解微分方程组，然后用非线性最小二乘模拟参数，急求回复，谢谢大家！ |
回复此楼» 猜你喜欢
B区学生调剂-兰州交通大学材料科学与工程学院
已经有12人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有13人回复
-
物理学I论文润色/翻译怎么收费?
已经有285人回复
-
桂林理工大学物理学专业招收调剂
已经有18人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生一名
已经有0人回复
-
【原创讨论】从电子约束到物质编辑：一套可迭代的环形磁场科技树
已经有0人回复
-
【方案分享】单环磁场+轴心控制+偏转导出电子束约束系统（可行性实验）
已经有6人回复
-
【修正版】单环用磁约束低速电子实验方案（简化版）
已经有0人回复
-
桂林理工大学物理学专业招收调剂，还有三个名额！！！
已经有20人回复
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
5