| 查看: 1446 | 回复: 7 | ||||
| 当前主题已经存档。 | ||||
liuyong2007金虫 (小有名气)
|
[交流]
【求助】用Matlab求解动力学方程参数的估计问题,程序运行有错误【已完成】
|
|||
|
哪位帮帮我,谢谢了.我用龙阁库塔积分和非线性最小二乘估计动力学参数,但是有错误.请高手指点一下,谢谢.程序如下: function KineticsEst1_int % 动力学方程为rA=dCA/dt=k1*CA*CB-k2*CC clear all clc global CAm CBm CCm t = [0 20 40 60 120 180 300]; CAm = [10 8 6 5 3 2 1]'; CBm = [8 7 6 5 3 2 1]'; CCm = [0 5 7 8 10 12 14]'; % 非线性拟合 beta0 = [0.0053 1.39]; tspan = [0 20 40 60 120 180 300]; CA0 = 10; [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,CA0,CAm,CBm,CCm) ci = nlparci(beta,resid,jacobian) % 拟合效果图(实验与拟合的比较) [t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1) tspan(end)],CA0,[],beta); plot(t,CAm,'bo',t4plot,CA4plot,'k-') legend('Exp','Model') xlabel('时间t, s') ylabel('浓度C_A, mol/L') % 参数辨识结果 fprintf('Estimated Parameters:\n') fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tk2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2)) % ------------------------------------------------------------------ function f = OptObjFunc(beta,tspan,CA0,CAm) [t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta); f = CAc - CAm; % ------------------------------------------------------------------ function dCAdt = KineticsEqs(t,CA,CB,CC,beta) dCAdt = beta(1)*CA*CB-beta(2)*CC; % k1= beta(1), k2= beta(2) [ Last edited by nono2009 on 2009-9-23 at 21:15 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 matlab |
» 猜你喜欢
有多少人是今天查系统知道结果的？
已经有19人回复
-
国社科又开始会评了，不知道这次命运如何
已经有9人回复
-
能否申诉？
已经有5人回复
-
为什么到现在没收到通知？
已经有5人回复
-
有没有仍没收到信息的
已经有8人回复
-
面上意见出来了
已经有18人回复
-
找导师
已经有7人回复
-
29号明天会评吗
已经有3人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有11人回复
-
我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？
已经有11人回复
allenhero1228
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1617.5
- 帖子: 201
- 在线: 29.7小时
- 虫号: 420486
- 注册: 2007-07-10
- 性别: GG
- 专业: 化工/二甲醚催化剂研究
2楼2008-12-11 12:33:19
liuyong2007
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 840
- 红花: 1
- 帖子: 169
- 在线: 82.4小时
- 虫号: 400205
- 注册: 2007-06-12
- 性别: GG
- 专业: 化学反应工程
mat;ab提示错误信息,希望高手给我改改
|
??? Error using ==> optim\private\lsqncommon User supplied function ==> OptObjFunc failed with the following error: Error using ==> KineticsEst1_int22>OptObjFunc Too many input arguments. Error in ==> lsqnonlin at 163 [x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ... Error in ==> KineticsEst1_int22 at 17 [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... |
3楼2008-12-11 12:41:34
allenhero1228
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1617.5
- 帖子: 201
- 在线: 29.7小时
- 虫号: 420486
- 注册: 2007-07-10
- 性别: GG
- 专业: 化工/二甲醚催化剂研究
4楼2008-12-11 13:06:55
snipher950
木虫 (正式写手)
- 应助: 2 (幼儿园)
- 金币: 3024.6
- 红花: 3
- 帖子: 797
- 在线: 22.5小时
- 虫号: 522854
- 注册: 2008-03-12
- 性别: MM
- 专业: 化工系统工程
5楼2008-12-11 14:20:27
snipher950
木虫 (正式写手)
- 应助: 2 (幼儿园)
- 金币: 3024.6
- 红花: 3
- 帖子: 797
- 在线: 22.5小时
- 虫号: 522854
- 注册: 2008-03-12
- 性别: MM
- 专业: 化工系统工程
★ ★ ★ ★ ★
woshilsh(金币+5,VIP+0):感谢回答！支持热心虫友，常来！呵呵！
woshilsh(金币+5,VIP+0):感谢回答！支持热心虫友，常来！呵呵！
|
可以运行了，改了很多地方，你自己琢磨一下。只不过结果不对，那是因为你的实验数据是虚假的，Matlab拟和不出来好的结果。 function KineticsEst1_int % 动力学方程为rA=dCA/dt=k1*CA*CB-k2*CC clear all clc global CAm CBm CCm t = [0 20 40 60 120 180 300]; CAm = [10 8 6 5 3 2 1]'; CBm = [8 7 6 5 3 2 1]'; CCm = [0 5 7 8 10 12 14]'; % 非线性拟合 beta0 = [0.0000053 10]; tspan = [0 20 40 60 120 180 300]; CA0 = 10; CB0 = 8; CC0 = 0; [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc,beta0,[0 0],[],[],tspan,[CA0 CB0 CC0]) ci = nlparci(beta,resid,jacobian) % 拟合效果图(实验与拟合的比较) [t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1) tspan(end)],[CA0 CB0 CC0],[],beta); plot(tspan,CAm,'bo',t4plot,CA4plot(:,1),'k-') legend('Exp','Model') xlabel('时间t, s') ylabel('浓度C_A, mol/L') % 参数辨识结果 fprintf('Estimated Parameters:\n') fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tk2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2)) % ------------------------------------------------------------------ function f = OptObjFunc(beta,tspan,CA0) global CAm CBm CCm [t Y] = ode45(@KineticsEqs,tspan,CA0,[],beta); if length(t)==length(tspan) f = Y - [CAm CBm CCm]; else f=100000; end % ------------------------------------------------------------------ function dy = KineticsEqs(t,Y,beta) CA=Y(1); CB=Y(2); CC=Y(3); dCAdt = beta(1)*CA*CB-beta(2)*CC; % k1= beta(1), k2= beta(2) dCBdt = beta(1)*CA*CB-beta(2)*CC; dCCdt = beta(2)*CC; dy=[dCAdt;dCBdt;dCCdt]; |
6楼2008-12-11 14:35:07
liuyong2007
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 840
- 红花: 1
- 帖子: 169
- 在线: 82.4小时
- 虫号: 400205
- 注册: 2007-06-12
- 性别: GG
- 专业: 化学反应工程
7楼2008-12-11 18:50:47
apolloking
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1086.1
- 散金: 35
- 红花: 2
- 帖子: 156
- 在线: 60.5小时
- 虫号: 75702
- 注册: 2005-06-20
- 专业: 能源化工
8楼2009-07-17 18:58:49
10