应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 未完成的matlab
51/1返回列表
查看: 791  |  回复: 5
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

805592836_01

铁虫 (初入文坛)

[求助] 未完成的matlab已有2人参与

有实验测得的一系列时间点
t=（0,180,300,420,600,900,1200,1800,2400）
CA=（1.845，1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124）

微分方程：k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA+k(3)^(1/6)*CB.^(1/2)).^7
CB=0.4为定值
希望由此得到三个参数k(1),k(2),k(3)，以及实验值与拟合值的比较图。现已经编了部分代码，但不是很确定对否，希望高手可以修改一下

clear all; clc
k0=[1,1,1];
lb=[0,0,0];
ub=[+inf,+inf,+inf];

tspan=[0,180,300,420,600,900,1200,1800,2400];

yexp=[1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124]';
y0=1.845;

% 使用fmincon()进行参数估计
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,lb,ub,option,tspan,y0,yexp);
fprintf('\tk1=%.4f\n',k(1)),  %本征动力学参数
fprintf('\tk2=%.4f\n',k(2)),  % A的吸附平衡常数KA
fprintf('\tk3=%.4f\n',k(3)),  % H2的吸附平衡常数KB
fprintf('The sum of the squares is:%.1e\n\n',fval),
k_fmincon=k;

% 使用函数Isqnonlin()进行参数估计
yy=[y0 yexp'];
options = optimset('MaxFunEvals',100000)
[k,resmorm,residual,exitflag,output,lambda,jacobian]=...
   Isqnonlin(@ObjFunc4LNL,k0,lb,ub,options,tspan,y0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('n\n使用Iaqunonlin()估计得到的参数值:\n'),Output

% 以函数fmincon()估计得到的结果为初值，使用Isqnonlin()进行参数估计
k0=k_fmincon;
[k,resmorm,residual,exitflag,output,lambda,jacobian]=...
   Isqnonlin(@ObjFunc4LNL,k0,lb,ub,options,tspan,y0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('n\n以fmincon()结果为初值，使用函数Isqnonlin()估计得到的参数值为:\n')
Output




%------------------------------------------------------------

function f=ObjFunc4Fmincon(k,tspan,y0,yexp)
[t y]=ode45(@KineticsEqs,tspan,y0,yexp,k);

%----------------------------------------------------------

function f=Objfunc4LNL(k,tspan,yexp)
[t y]=ode45(@KineticEqs,tspan,tspan,[],k);



%----------------------------------------------------------
function dydt = KineticsEqs(t,y,k)
y0=1.845;
CB=y;
CB0=0.4;
CA0=y0;
k1=k(1);
KA=k(2);
KB=k(3);
dxdt=k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA0+k(3)^(1/6)*CB.^(1/2)).^7
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
加油加油加油……
1楼2014-03-12 08:57:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

805592836_01

铁虫 (初入文坛)
引用回帖:
4楼: Originally posted by dingd at 2014-03-12 10:35:53
不可能有好结果，检查下模型吧。...

嗯，好的吧！同样感谢
赞一下 回复此楼
加油加油加油……
5楼2014-03-12 11:31:29
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
k1, k2, k3要求都大于0的话效果会很差的。
赞一下 回复此楼
2楼2014-03-12 09:56:13
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

805592836_01

铁虫 (初入文坛)
引用回帖:
2楼: Originally posted by dingd at 2014-03-12 09:56:13
k1, k2, k3要求都大于0的话效果会很差的。

嗯，都要大于0的……
赞一下 回复此楼
加油加油加油……
3楼2014-03-12 10:02:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)
引用回帖:
3楼: Originally posted by 805592836_01 at 2014-03-12 10:02:43
嗯，都要大于0的……...

不可能有好结果，检查下模型吧。
赞一下 回复此楼
4楼2014-03-12 10:35:53
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭