| 查看: 292 | 回复: 0 | ||
Shirley89银虫 (小有名气)
|
[求助]
求大神 帮忙，帮我看看问题出在哪里？
|
|
function parafit % % r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c)); % % dCAdt = - r1; clear all clc % t/min CA / mol/L Kinetics=[0 7.029559717 1 3.960706476 2 3.163862085 3 2.963057756 6 2.820434159 9 2.698109429 12 2.556212898 15 2.447100967 18 2.207246045 21 2.13921916 31 2.000403189 41 1.917172016 51 1.851299168 61 1.636287388]; c=Kinetics( :,1)'; k0 = [0.671 0.02253 0.641 0.04028]; % 参数初值 lb = [0 0 0 0]; % 参数下限 ub = [100 100 100 100]; % 参数上限 x0 = [7.029559717]; yexp = Kinetics; warning off % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.11f\n',k(1)) fprintf('\tk2 = %.11f\n',k(2)) fprintf('\tk3 = %.11f\n',k(3)) fprintf('\tk4 = %.11f\n',k(4)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fm= k; warning off % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.11f\n',k(1)) fprintf('\tk2 = %.11f\n',k(2)) fprintf('\tk3 = %.11f\n',k(3)) fprintf('\tk4 = %.11f\n',k(4)) fprintf(' The sum of the squares is: %.1e\n\n',resnorm) k_ls = k; output warning off % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fm; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.11f\n',k(1)) fprintf('\tk2 = %.11f\n',k(2)) fprintf('\tk3 = %.11f\n',k(3)) fprintf('\tk4 = %.11f\n',k(4)) fprintf(' The sum of the squares is: %.1e\n\n',resnorm) k_fmls = k; output tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61]; [t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls); figure; plot(t,x(:,1),t,yexp(:,2),'*');legend('ca-pr','ca-real') x(:,1) figure;plot(t,x); hold on % ------------------------------------------------------------------ function f = ObjFunc7Fmincon(k,x0,yexp) tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,2) = x(:,1); f = sum((y(:,2)-yexp(:,2)).^2) ; % ------------------------------------------------------------------ function f = ObjFunc7LNL(k,x0,yexp) tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,2) = x(:,1); f1 = y(:,2) - yexp(:,2); f = [f1]; % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) Kinetics=[0 7.029559717 1 3.960706476 2 3.163862085 3 2.963057756 6 2.820434159 9 2.698109429 12 2.556212898 15 2.447100967 18 2.207246045 21 2.13921916 31 2.000403189 41 1.917172016 51 1.851299168 61 1.636287388]; c=Kinetics(:,1)'; for i=1:14 c(i)=c(:,i); r1 =(k(1)*k(2)*c(i)^2-0.25*k(3)*k(4)*(9.404-c(i))^2)/(1+k(2)*c(i)+0.5*k(4)*(9.404-c(i))); dCAdt = - r1; dxdt = [dCAdt]; end 运行出来是图中那样，不知道问题出在哪里，改变初值也不行。  H8_DJFZ~6NZC}07J0EXL1]L.png |
回复此楼» 猜你喜欢
如何从铁电相到顺电相。
已经有1人回复
-
各位大佬，求一份最新Wien2k版本作为学习用途
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有248人回复
-
南科大活性流体和软物质课题组诚招2027级博士、硕士研究生和博后
已经有42人回复
-
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有24人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
找到一些相关的精华帖子，希望有用哦~
求大神 帮忙解决加氢的问题 老是出问题 附件在下方
已经有7人回复
- 求助打算帮我看看我的引物二聚体很亮，目的条带很暗怎么破
已经有6人回复
- 求大神给我看看这个反应咋没反应
已经有4人回复
- 求大神帮忙拟合一个非线性方程，求出模型参数
已经有15人回复
- 关于UDS的疑问
已经有8人回复
- 很简单的Fortran赋值出错啦，求大神们帮忙找bug...
已经有11人回复
- 按照ansys官方教程一步步练习，求解时经常出现问题，求大神分析
已经有14人回复
- 帮我看看问题出在哪？
已经有12人回复
- 英国商务签证 学术会议 遭到拒签 求好心人帮忙看看问题出在哪里
已经有23人回复
- 求大神帮忙看一下断口
已经有6人回复
- 审稿意见回来了 有些理解不是很透彻啊 求各位大神帮忙看看是什么意思啊？？？
已经有20人回复
- 哪位好心人看我哪里出问题了
已经有4人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
100