24小时热门版块排行榜    

查看: 657  |  回复: 3

zxqworld

新虫 (初入文坛)

[求助] 求大神帮忙看一下 我写的积分法 求取动力学参数的 程式 已有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个,还请各位大神帮忙看一下。。万分感谢!!
回复此楼

» 猜你喜欢

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

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

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
微分方程拟合你试试1stOpt,更简单好用些。
2楼2014-05-23 10:21:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zxqworld

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by dingd at 2014-05-23 10:21:25
微分方程拟合你试试1stOpt,更简单好用些。

好的 ,我试试看 。因为比较急 ,没用过 ,不知道会不会比较麻烦。感谢!
3楼2014-05-23 10:43:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

云少爷

新虫 (初入文坛)

【答案】应助回帖

楼主可以看一下实用化工计算机模拟-Matlab在化学工程中的应用,第七章讲到了你问的问题,我记得有一个跟你问题类似的程序
4楼2016-11-09 11:58:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zxqworld 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 化学工程321分求调剂 +18 大米饭! 2026-03-15 22/1100 2026-03-21 20:20 by HH领袖
[考研] 本人考085602 化学工程 专硕 +20 不知道叫什么! 2026-03-15 22/1100 2026-03-21 19:03 by ColorlessPI
[考研] 求调剂 +3 13341 2026-03-20 3/150 2026-03-21 18:28 by 学员8dgXkO
[考研] 【考研调剂】化学专业 281分,一志愿四川大学,诚心求调剂 +11 吃吃吃才有意义 2026-03-19 11/550 2026-03-21 18:23 by 学员8dgXkO
[考研] 求助 +5 梦里的无言 2026-03-21 6/300 2026-03-21 17:51 by 学员8dgXkO
[考研] 一志愿西安交通大学材料工程专业 282分求调剂 +8 枫桥ZL 2026-03-18 10/500 2026-03-21 15:29 by Shawn0911
[考研] 材料与化工(0856)304求 B区 调剂 +3 邱gl 2026-03-21 3/150 2026-03-21 13:47 by lature00
[考研] 070300化学319求调剂 +7 锦鲤0909 2026-03-17 7/350 2026-03-21 03:46 by JourneyLucky
[考研] 一志愿天津大学化学工艺专业(081702)315分求调剂 +12 yangfz 2026-03-17 12/600 2026-03-21 03:30 by JourneyLucky
[考研] 华东师范大学-071000生物学-293分-求调剂 +3 研究生何瑶明 2026-03-18 3/150 2026-03-21 01:30 by JourneyLucky
[考研] 22408 344分 求调剂 一志愿 华电计算机技术 +4 solanXXX 2026-03-20 4/200 2026-03-20 23:49 by alg094825
[考研] 一志愿西南交通 专硕 材料355 本科双非 求调剂 +5 西南交通专材355 2026-03-19 5/250 2026-03-20 21:10 by JourneyLucky
[考研] 一志愿吉林大学材料学硕321求调剂 +11 Ymlll 2026-03-18 15/750 2026-03-20 19:40 by 丁丁*
[考研] 261求B区调剂,科研经历丰富 +3 牛奶很忙 2026-03-20 4/200 2026-03-20 19:34 by JourneyLucky
[考研] 环境工程调剂 +9 大可digkids 2026-03-16 9/450 2026-03-20 17:38 by 醉在风里
[考研] 081700化工学硕调剂 +3 【1】 2026-03-16 3/150 2026-03-19 23:40 by edmund7
[考研] 材料,纺织,生物(0856、0710),化学招生啦 +3 Eember. 2026-03-17 9/450 2026-03-18 10:28 by Eember.
[考研] 东南大学364求调剂 +5 JasonYuiui 2026-03-15 5/250 2026-03-16 21:28 by 木瓜膏
[考研] 333求调剂 +3 文思客 2026-03-16 7/350 2026-03-16 18:21 by 文思客
[考研] 304求调剂 +5 素年祭语 2026-03-15 5/250 2026-03-16 17:00 by 我的船我的海
信息提示
请填处理意见