24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 3631  |  回复: 15

lxyy

银虫 (小有名气)

[求助] matlab拟合方程参数时初值的选择

A+B=C  A+W=R  R+B=C+W
先考虑的三个反应都是可逆反应,有以下几个式子:编写程序如下
%30℃
clear all
clc
global keq1 keq2 t0 c0

keq1 = 3.94;         % 平衡常数k1
keq3 = 1.08;   % 平衡常数k3

t0 = [5,10, 15, 30, 45, 60,90,120,160,200,240,300,360];
c0 = [5.939601,5.203973,0,0.022484,0.078956]        % A B C R W 组成的初始浓度
ci=[3.742508811        4.213566482        0.990406938        0.252498405        0.090415512;
3.045716838        3.387013322        1.816960097        0.215546741        0.081072022;
2.504754523        2.901991595        2.301981824        0.193427102        0.099462604;
1.65098702        2.036686653        3.167286766        0.170236582        0.095144044;
1.271633833        1.655460403        3.548513016        0.124730781        0.10334349;
1.129887499        1.630815394        3.573158025        0.129798695        0.097764543;
1.01239208        1.460748979        3.74322444        0.144649391        0.073614958;
1.007417603        1.429180292        3.774793127        0.123477124        0.094249307;
0.990005723        1.499216939        3.704756481        0.135524067        0.076072022;
0.970073726        1.493966669        3.71000675        0.130791851        0.076387812;
0.96903024        1.429705429        3.77426799        0.119510017        0.101202216;
0.980873385        1.457988124        3.745985296        0.124189274        0.097022161;
0.98477262        1.453402965        3.750570455        0.1200936         0.0762608];

% cA cC cR 对应时间数值

k0=[1 1 1];   % k1+, k2, k3+ 初始向值
lb = [0 0 0];
ub = [1000  1000 500];  % 上下限

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag] = lsqnonlin(@ObjFunc,k0,lb,ub,[],ci);      

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

% ------------------------------------------------------------------
function f = ObjFunc(k,ci)    % 目标函数
global t0  c0

[t,c_cal] = ode45(@Euqations,t0,c0,[],k);
f= sum((c_cal-ci)^2);

% ------------------------------------------------------------------
function dcdt = Euqations(t,c,k)
global keq1 keq2

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

cA=c(:,1);
cB=c(:,2);
cC=c(:,3);
cR=c(:,4);
cW=c(:,5);

dcAdt = k1plus*cA*cB-k1minus*cC+k2*cA*cW;
dcCdt = -k1plus*cA*cB+k1minus*cC+k3plus*cB*cR-k3minus*cC*cW;
dcRdt = k2*cA*cW-k3plus*cB*cR+k3minus*cC*cW;
dcdt = [dcAdt; dcCdt;dcRdt];
但是k的初值一直影响结果,程序要怎么改一下初值对结果影响不大(这个程序有问题,)
回复此楼

» 猜你喜欢

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

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

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
微分方程dcBdt, dcWdt的表达式呢?
2楼2012-03-28 10:37:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
2楼: Originally posted by dingd at 2012-03-28 10:37:30:
微分方程dcBdt, dcWdt的表达式呢?

我写在附件里面了

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 太残忍doc.doc
  • 2012-03-28 10:52:17, 30.5 K
3楼2012-03-28 10:52:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖


zhangguangping: 金币+1, 谢谢指导! 2012-03-28 14:16:22
最好把你的专业问题用通俗的数学语言来描述,非化工专业的很难理解的;
你的问题应该就是一微分方程拟合问题,仅需描述清楚:

1:微分方程组表达式
2:对应的数据
3:要求的参数
4楼2012-03-28 11:07:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
4楼: Originally posted by dingd at 2012-03-28 11:07:18:
最好把你的专业问题用通俗的数学语言来描述,非化工专业的很难理解的;
你的问题应该就是一微分方程拟合问题,仅需描述清楚:

1:微分方程组表达式
2:对应的数据
3:要求的参数

不好意思,我的微分方程组就是附件里面标(1)(2)(3)的三个式子,对应的数据为

时间:t0 = [5,10, 15, 30, 45, 60,90,120,160,200,240,300,360];
不同时刻对应的数据
c0 = [5.939601,5.203973,0,0.022484,0.078956]        % A B C R W 组成的初始浓度
ci=[3.742508811        4.213566482        0.990406938        0.252498405        0.090415512;
3.045716838        3.387013322        1.816960097        0.215546741        0.081072022;
2.504754523        2.901991595        2.301981824        0.193427102        0.099462604;
1.65098702        2.036686653        3.167286766        0.170236582        0.095144044;
1.271633833        1.655460403        3.548513016        0.124730781        0.10334349;
1.129887499        1.630815394        3.573158025        0.129798695        0.097764543;
1.01239208        1.460748979        3.74322444        0.144649391        0.073614958;
1.007417603        1.429180292        3.774793127        0.123477124        0.094249307;
0.990005723        1.499216939        3.704756481        0.135524067        0.076072022;
0.970073726        1.493966669        3.71000675        0.130791851        0.076387812;
0.96903024        1.429705429        3.77426799        0.119510017        0.101202216;
0.980873385        1.457988124        3.745985296        0.124189274        0.097022161;
0.98477262        1.453402965        3.750570455        0.1200936         0.0762608];
要求的是ki+,k2,k3+,
另外方程中的ki-=ki+/3.14,k3-=k3+/1.08

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 太残忍doc (1)1.doc
  • 2012-03-28 12:18:21, 26 K
5楼2012-03-28 12:18:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

数据项B、W有对应的微分方程吗?还是自由变量?
6楼2012-03-28 14:58:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
6楼: Originally posted by dingd at 2012-03-28 14:58:43:
数据项B、W有对应的微分方程吗?还是自由变量?

有的,我写在附件里了,真是太谢谢了

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 太残忍doc (1)1.doc
  • 2012-03-28 15:05:19, 32.5 K
7楼2012-03-28 15:05:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

怎么感觉和挤牙膏似的!!

公式写的对吗?是否好好检查过?和1楼的代码有不少不同,比如对r2的定义,"dCB/wdt"应该是“dCB/dt”吗?公式(2)和(5)都给的是“dCB/dt”,那“dCC/dt”呢?
8楼2012-03-28 16:30:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
8楼: Originally posted by dingd at 2012-03-28 16:30:05:
怎么感觉和挤牙膏似的!!

公式写的对吗?是否好好检查过?和1楼的代码有不少不同,比如对r2的定义,"dCB/wdt"应该是“dCB/dt”吗?公式(2)和(5)都给的是“dCB/dt”,那“dCC/dt”呢?

dCB/wdt"应该是“dCB/dt,没错
对不起我太心急了,公式已经检查过,我保证这次没错了,真的

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 太残忍doc(1)2.doc
  • 2012-03-28 16:38:07, 34.5 K
9楼2012-03-28 16:38:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

供参考:
CODE:
Constant keq1 = 3.94, keq3 = 1.08;
ConstStr k1minus= k1plus/keq1,k3minus= k3plus/keq3,
         r1=k1plus*cA*Cb-k1minus*cC,
         r2=k2*cA*cW,
         r3=k3plus*cB*cR-k3minus*cC*cW;
Variable t,ca,cb,cc,cr,cw;
ODEFunction cA' = -r1-r2;
            cB' = -r1-r3;
            cC' = r1+r3;
            cR' = r2-r3;
            cW' = r3-r2;
Data;
//t,cA,cB,cC,cR,cW
0        5.939601        5.203973        0        0.0282484        0.078956
5        3.742508811        4.213566482        0.990406938        0.252498405        0.090415512
10        3.045716838        3.387013322        1.816960097        0.215546741        0.081072022
15        2.504754523        2.901991595        2.301981824        0.193427102        0.099462604
30        1.65098702        2.036686653        3.167286766        0.170236582        0.095144044
45        1.271633833        1.655460403        3.548513016        0.124730781        0.10334349
60        1.129887499        1.630815394        3.573158025        0.129798695        0.097764543
90        1.01239208        1.460748979        3.74322444        0.144649391        0.073614958
120        1.007417603        1.429180292        3.774793127        0.123477124        0.094249307
160        0.990005723        1.499216939        3.704756481        0.135524067        0.076072022
200        0.970073726        1.493966669        3.71000675        0.130791851        0.076387812
240        0.96903024        1.429705429        3.77426799        0.119510017        0.101202216
300        0.980873385        1.457988124        3.745985296        0.124189274        0.097022161
360        0.98477262        1.453402965        3.750570455        0.1200936        0.0762608

均方差(RMSE): 0.384192283597719
残差平方和(SSE): 9.59424120044197
相关系数(R): 0.998354081222552
相关系数之平方(R^2): 0.996710871493726
决定系数(DC): 0.995560449063788

参数                  最佳估算
--------------------        -------------
k1plus        0.0122448530178222
k2        0.044692817673755
k3plus        -0.0135416873776772


10楼2012-03-28 17:10:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 lxyy 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 375求调剂 +6 雨夏整夜 2026-03-29 6/300 2026-03-30 10:21 by herarysara
[考研] 一志愿北京化工大学材料与化工(085600)296求调剂 +17 稻妻小编 2026-03-26 17/850 2026-03-29 19:07 by 学员8dgXkO
[考研] 289求调剂 +5 BrightLL 2026-03-29 5/250 2026-03-29 17:24 by zhyzzh
[考研] 一志愿南昌大学324求调剂 +5 hanamiko 2026-03-29 5/250 2026-03-29 16:46 by 学员8dgXkO
[考研] 070305高分子化学与物理 304分求调剂 +12 c297914 2026-03-28 12/600 2026-03-29 16:21 by Serene1974
[考研] 085602 化工专硕 338分 求调剂 +12 路痴小琪 2026-03-27 12/600 2026-03-28 15:41 by L135790
[考研] 一志愿厦门大学化学学硕307求调剂 +10 y7czhao 2026-03-26 10/500 2026-03-28 14:23 by 唐沐儿
[考研] 0856,材料与化工321分求调剂 +12 大馋小子 2026-03-27 13/650 2026-03-28 10:56 by self2008
[考研] 339求调剂,想调回江苏 +6 烤麦芽 2026-03-27 8/400 2026-03-28 10:40 by 烤麦芽
[考研] 求调剂推荐 材料 304 +15 荷包蛋hyj 2026-03-26 15/750 2026-03-28 04:13 by fmesaito
[考研] 265求调剂11408 +3 刘小鹿lu 2026-03-27 3/150 2026-03-27 20:53 by nihaoar
[考研] 化学调剂 +4 爱吃番茄的旭 2026-03-24 5/250 2026-03-27 17:50 by kiokin
[考研] 287求调剂 +10 land xuxu 2026-03-26 10/500 2026-03-27 15:33 by 帕尔马拉特
[考研] 279 分 求调剂 +4 睡个好觉_16 2026-03-24 4/200 2026-03-27 15:05 by 醉在风里
[考研] 08开头275求调剂 +4 拉谁不重要 2026-03-26 4/200 2026-03-27 14:12 by Delta2012
[考研] 304材料求调剂 +4 钟llll 2026-03-26 4/200 2026-03-27 03:42 by wxiongid
[考研] 351求调剂 +4 麦克阿磊 2026-03-24 4/200 2026-03-27 00:32 by wxiongid
[考研] 打过很多竞赛,085406控制工程300分,求调剂 +3 askeladz 2026-03-26 3/150 2026-03-26 09:08 by 给你你注意休息
[考研] 334分 一志愿武理-080500 材料求调剂 +4 李李不服输 2026-03-25 4/200 2026-03-25 21:26 by 星空星月
[考研] 一志愿吉大化学322求调剂 +4 17501029541 2026-03-23 6/300 2026-03-24 10:21 by 戴围脖的小蚊子
信息提示
请填处理意见