24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

805592836_01

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 823
帖子: 10
在线: 23.4小时
虫号: 2393210
注册: 2013-04-01
性别: MM
专业: 能源化工

[求助] 未完成的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人回复
求取一些关于纳米材料和纳米技术相关的英文PPT。已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有83人回复
【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有19人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
26申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复

» 本主题相关价值贴推荐，对您同样有帮助:

matlab中的计算问题？已经有8人回复
MATLAB求助，关于importdata的问题已经有6人回复
用matlab制作的M文件保存在哪里才能直接调用？已经有9人回复
matlab一运行就出现的问题已经有7人回复
matlab的m文件关联问题已经有17人回复
求助一个matlab关于分类讨论的问题已经有9人回复
matlab编程中遇到的问题已经有4人回复
帮忙修改Matlab中出现的错误已经有4人回复
matlab读取图像的问题已经有12人回复
分享matlab的一些代码（入门贴）已经有19人回复
求MATLAB清华版的教学视频已经有14人回复
关于Matlab作图的问题已经有14人回复
MATLAB初学者，这样的式子怎么就是画不了图，总是出错!求指导!!!!!!!!!!!!! 已经有8人回复
matlab程序将图中的数据记录下来已经有4人回复
Matlab在地学中的应用已经有24人回复
matlab中dydt的编写要求已经有12人回复
求助svm在matlab中的应用已经有7人回复
关于matlab 2010b 的一个小问题已经有7人回复
MATLAB安装后出现一堆java的问题已经有12人回复
【求助】关于matlab的高维数组已经有6人回复

加油加油加油……

1楼2014-03-12 08:57:58

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

专家经验: +1059
计算强帖: 8
应助: 1712 (讲师)
贵宾: 8.888
金币: 68120.7
散金: 1938
红花: 443
沙发: 4
帖子: 4373
在线: 3291.4小时
虫号: 1122189
注册: 2010-10-14
专业: 宇宙学
管辖: 计算模拟区

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
805592836_01: 金币+15, ★★★很有帮助, 看来还是模型不合适，三个参数要都大于0，不过还是感谢大神的帮忙。 2014-03-12 21:06:07

原代码问题很多，比如MATLAB中不存在Isqnonlin函数，系拼写错误，正确的为 lsqnonlin；此外用fmincon函数算出来的结果输入给lsqnonlin函数作初值意义不是很大，单独调用lsqnonlin即可。
根据实验数据可知，数据呈单调递减的趋势，拟合公式是关于时间的导数，应小于0，所以k1 k2 k3不应该全部大于0，否则不可能有号的拟合效果。
在以下给出程序中，把k1的取值范围可为0，具体代码以及计算结果如下。

CODE:

function test_123

clear all;clc

format long

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;

k0=[-1 1 1];

lb=-[inf 0 0];

ub=+[inf inf inf];        

yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

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',resnorm)

ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);

[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k); 

y=XXsim(2:end);

xexp=yexp;

R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);

fprintf('\n\t相关系数之平方R^2 = %.6f',R2);

plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best'),

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

function f = ObjFunc(k,tspan,y0,yexp)           % 目标函数

[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;

ysim = Xsim(2:end);

f=ysim-yexp;

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

function dxdt = 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;