| 查看: 867 | 回复: 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求解微分方程组，然后用非线性最小二乘模拟参数，急求回复，谢谢大家！ |
回复此楼» 猜你喜欢
浙大物院表面物理课题组长期招收普博及直博学生
已经有259人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有47人回复
-
物理学I论文润色/翻译怎么收费?
已经有137人回复
-
25届硕士申请，材料计算方向
已经有2人回复
-
课题组招1-2名硕士
已经有1人回复
-
机时越来越贵之后，我把一部分 DFT 流程搬到了 AI 计算平台上
已经有0人回复
-
斯坦福大学与Elsevier世界前2%学者新排名查询
已经有1人回复
-
27申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有19人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
-
apl投稿状态
已经有2人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有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
50