| 查看: 1195 | 回复: 4 | ||
| 【悬赏金币】回答本帖问题，作者Lulahey将赠送您 5 个金币 | ||
[求助]
请会MATLAB的大佬帮我看一下动力学模型存在的问题已有1人参与
|
||
|
我看了一篇和我的体系很像的文献，里面动力学用了LHHW的动力学模型，我就按照他的模型写了MATLAB 的代码 function nihe_0316 % 微分方程的参数估计 % clear all; clc; format long; global t0; % k10,E,Ka 对应k(1)、k(2)、k(3) k0 =[5.33E9;93.06;0.0437]; % 参数初值，重要！，可将结果反复代入迭代计算！！ x0 = 0.63713; %初始状态，没有就取第一行数据 lb = [0 0 0]; % 参数下限,根据拟合结果，适当调整 ub = [+inf +inf +inf]; % 参数上限 % 实验数据: expdata= ... [0 0.63713 60 0.79642 120 0.87377 240 0.90782 360 0.91735]; t0=expdata(:,1); yexp =expdata(:,2); % yexp: 对应实验数据[Xb] % ------------------------------------------------------------------ % 使用函数lsqnonlin()进行参数估计 k=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk10 = %.8f\n',k(1)) fprintf('\tE = %.8f\n',k(2)) fprintf('\tKa = %.8f\n',k(3)) % ---------------------------------------------------- tspan = [t0(1) t0(end)]; [t,xx] = ode45(@KineticEqs,tspan,x0,[],k); x=spline(t,xx,t0); %插值，使理论值和实验值个数一致 % 计算相关系数 A=corrcoef(yexp,x(:,1)); fprintf('\nXb的相关系数R为:\n') fprintf('\tR = %.8f\n',A(1,2)) % % 计算决定系数 r(1)=1-sum((yexp-x(:,1)).^2)./sum((yexp-mean(yexp)).^2); % ********************************决定系数 fprintf('\nXb的决定系数R^2为:\n') fprintf('\tR^2 = %.8f\n',r(1)) % 计算均方根误差 rmse(1)=sqrt(sum((yexp-x(:,1)).^2)/length(t0)); fprintf('\nXb的均方根误差RMSE为:\n') fprintf('\tRMSE = %.8f\n',rmse(1)) % % 计算残差平方和 SSE（Sum of Squares for Error） ssr(1)=sum((yexp(:,1)-x(:,1)).^2); fprintf('\nXb的残差平方和为:\n') fprintf('\tSSE = %.8f\n',ssr(1)) % ---------------------------------------------------- figure(1) plot(t,xx(:,1),'b-',t0,yexp,'ro'),legend('理论值','实验值','Location','best') xlabel('时间');ylabel('实验值N'); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = [t0(1) t0(end)]; [t,xx] = ode45(@KineticEqs,tspan,x0,[],k); x=spline(t,xx,t0); f = x - yexp; % 理论值与实验值的差值，残差 end % ------------------------------------------------------------------ function dXbdt = KineticEqs(t,Xb,k) % 微分方程, k10,E,Ka 对应k(1)、k(2)、k(3) T=358.15; m=3; N0=0.1467; Cb0=0.5564; Ca=Cb0*(3-Xb); Cb=Cb0*(1-Xb); Cc=Cb0*Xb; Keq=exp(4873/T-13.813); Kb=5.2112*k(3); Kc=2.2695*k(3); dXbdt = (m*k(1)*exp(-k(2)/8.314*T)*(Ca*Cb-Cc/Keq))/(N0*(1+k(3)*Ca+Kb*Cb+Kc*Cc)^2); end % ------------------------------------------------------------------ end 最后X的计算值与t无关，书上有一个例子是dc/dt=f（c）用了RT+最小二乘法模拟了两个参数，我这个也是dx/dt=f（x）用RT+最小二乘法模拟三个参数就不对了，有高手可以指点一下吗？  QQ截图20200317002352.png |
回复此楼» 猜你喜欢
博士白读了
已经有53人回复
-
with editor 两个月了，什么原因？
已经有18人回复
-
优青已经开始了么？
已经有3人回复
-
请问UV灯是365nm的，那么选光引发剂的波长选多少的？要完全一致吗？
已经有5人回复
-
江苏理工学院要慎选
已经有3人回复
-
博士招生-211农业院校
已经有7人回复
-
材料方向24博士申请/一作SCI三篇
已经有5人回复
-
失之东隅 收之桑榆
已经有6人回复
-
真的好想读博！
已经有23人回复
-
基金开始函评了吗？
已经有17人回复
七月民谣
木虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1095.8
- 散金: 3088
- 红花: 3
- 帖子: 963
- 在线: 55.8小时
- 虫号: 5864427
- 注册: 2017-03-07
- 专业: 高分子材料的加工与成型
2楼2020-06-10 07:42:29
3楼2024-03-18 10:49:18
johnliu1983
至尊木虫 (著名写手)
- 应助: 33 (小学生)
- 金币: 13249.7
- 红花: 7
- 帖子: 1778
- 在线: 121.2小时
- 虫号: 370817
- 注册: 2007-05-14
- 性别: GG
- 专业: 细胞生物学
4楼2024-03-20 09:11:10
johnliu1983
至尊木虫 (著名写手)
- 应助: 33 (小学生)
- 金币: 13249.7
- 红花: 7
- 帖子: 1778
- 在线: 121.2小时
- 虫号: 370817
- 注册: 2007-05-14
- 性别: GG
- 专业: 细胞生物学
5楼2024-03-20 09:36:59
15