| 查看: 656 | 回复: 3 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
求大神帮忙看一下 我写的积分法 求取动力学参数的 程式 已有2人参与
|
||
|
function KineticsEst1_int % 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n clear all clc global X t =[30:30:120]; X1=[67.37 58.28 45.32 32.1]; X2=[7.88 17.88 17.45 26.43]; X3=[2.48 9.09 11.84 13.65]; X4=[1.58 0 0.94 0]; X5=[20.69 14.75 24.45 27.82]; % 非线性拟合 k0 = [0.00594 0.00391 0.00222 0.00810 0.00020 0.01071 ]; n0 =[1 1 1 1.43 1.38 1.26]; beta0=[k0; n0]; tspan = [30:120]; X0 = [67.37 7.88 2.48 1.58 20.69]; [beta,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,X0,X); ci = nlparci(beta,resid,jacobian); % 拟合效果图(实验与拟合的比较) figure(1) [t4plot X4plot] = ode45(@KineticsEqs,[tspan],X0,[],beta); plot(t,X(:,1),'bo',t4plot,X4plot,'k-') legend('Exp','Model') xlabel('时间 min') ylabel('收率C_C, %') % 残差关于拟合值的残差图 [t Xc] = ode45(@KineticsEqs,tspan,X0,[],k,n); figure plot(Xc,resid,'*') xlabel('收率拟合值（%）') ylabel('残差R (%)') refline(0,0) % 参数辨识结果 fprintf('Estimated Parameters:\n') fprintf('\tk = %.4f ± %.4f\n',k,ci(1,2)-k) fprintf('\tn = %.2f ± %.2f\n',n,ci(2,2)-n) % ------------------------------------------------------------------ function f = OptObjFunc(k,n,tspan,X0,X) tspan=[30:30:120] [t Xc] = ode45(@KineticsEqs,tspan,X0,[],k,n); f1 = Xc(:,1) - X(:,1); f2 = Xc(:,2) - X(:,2); f3 = Xc(:,3) - X(:,3); f4 = Xc(:,4) - X(:,4); f5 = Xc(:,5) - X(:,5); f = f1+ f2+ f3+ f4+ f5; % ------------------------------------------------------------------ function dXdt = KineticsEqs(t,X,k,n) dXdt = ... [ ( -(k(1)+k(2)+k(3)+k(4))*X(1)^n(1) ) ( k(2)*X(1)^n(2)) ( k(1)*X(1)^n(1)+k(5)*X(4)^n(5)+k(6)*X(5)^n(6) ) ( k(3)^X(1)^n(3)-k(5)*X(4)^n(5)) ( k(4)*X(1)^n(4)-k(6)*X(5)^n(6)) ]; 程序错误好多，但是不知道怎么改，尤其是 非线性拟合部分，关于矩阵beta，我是参照一本书上的一个例题写的，但是他只有一个方程，我有5个，还请各位大神帮忙看一下。。万分感谢！！ |
回复此楼» 猜你喜欢
336求调剂
已经有5人回复
-
277材料科学与工程080500求调剂
已经有6人回复
-
333求调剂
已经有4人回复
-
286求调剂
已经有8人回复
-
298求调剂
已经有4人回复
-
306求0703调剂一志愿华中师范
已经有5人回复
-
299求调剂
已经有7人回复
-
26考研一志愿中国石油大学(华东)305分求调剂
已经有6人回复
-
初试 317
已经有5人回复
-
296求调剂
已经有3人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 航空动力学报审稿求大神帮助
已经有3人回复
- 求ANSYS大神，热力学问题，我都快崩溃了
已经有8人回复
- 求助各位大神量子力学入门书籍，请知道的进来说下
已经有28人回复
4楼2016-11-09 11:58:17
dingd
铁杆木虫 (职业作家)
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2014-05-23 10:21:25
3楼2014-05-23 10:43:51
10