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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab拟合方程参数时初值的选择
161/2返回列表12下一页
查看: 3608  |  回复: 15
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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的初值一直影响结果，程序要怎么改一下初值对结果影响不大（这个程序有问题，）
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2012-03-28 10:24:19
已阅   回复此楼   关注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 的主题更新
161/2返回列表12下一页
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 南京大学化学376求调剂 +3 hisfailed 2026-03-19 6/300 2026-03-20 23:43 by hisfailed
[考研] 南昌大学材料专硕311分求调剂 +5 77chaselx 2026-03-20 5/250 2026-03-20 23:42 by lovewei0727
[考研] 一志愿西南交通 专硕 材料355 本科双非 求调剂 +5 西南交通专材355 2026-03-19 5/250 2026-03-20 21:10 by JourneyLucky
[考研] 295材料求调剂，一志愿武汉理工085601专硕 +5 Charlieyq 2026-03-19 5/250 2026-03-20 20:35 by JourneyLucky
[考研] 279分求调剂 一志愿211 +10 chaojifeixia 2026-03-19 11/550 2026-03-20 20:32 by JourneyLucky
[考研] 08工学调剂 +5 用户573181 2026-03-20 5/250 2026-03-20 15:47 by xia_2003
[考研] 281求调剂（0805） +14 烟汐忆海 2026-03-16 25/1250 2026-03-20 15:47 by yuncha
[考研] 304求调剂 +5 曼殊2266 2026-03-18 5/250 2026-03-20 09:00 by ZHANG0tao
[考研] 307求调剂 +9 冷笙123 2026-03-17 9/450 2026-03-19 22:44 by 学员8dgXkO
[考研] 梁成伟老师课题组欢迎你的加入 +9 一鸭鸭哟 2026-03-14 11/550 2026-03-19 17:22 by ！本暗一次！
[考研] 一志愿天津大学化学工艺专业（081702）315分求调剂 +11 yangfz 2026-03-17 11/550 2026-03-19 15:06 by houyaoxu
[考研] 一志愿武理材料305分求调剂 +5 想上岸的鲤鱼 2026-03-18 6/300 2026-03-18 17:53 by 无际的草原
[考研] 297求调剂 +8 戏精丹丹丹 2026-03-17 8/400 2026-03-18 14:30 by laoshidan
[考研] 278求调剂 +5 烟火先于春 2026-03-17 5/250 2026-03-18 08:43 by 星空星月
[考研] 334求调剂 +3 志存高远意在机� 2026-03-16 3/150 2026-03-18 08:34 by lm4875102
[考研] 一志愿，福州大学材料专硕339分求调剂 +3 木子momo青争 2026-03-15 3/150 2026-03-17 07:52 by laoshidan
[考研] 11408 一志愿西电，277分求调剂 +3 zhouzhen654 2026-03-16 3/150 2026-03-17 07:03 by laoshidan
[基金申请] 今年的国基金是打分制吗？ 50+3 zhanghaozhu 2026-03-14 3/150 2026-03-16 17:07 by 北京莱茵润色
[考研] 304求调剂 +4 ahbd 2026-03-14 4/200 2026-03-16 16:48 by 我的船我的海
[考研] 070305求调剂 +3 mlpqaz03 2026-03-14 4/200 2026-03-15 11:04 by peike
信息提示
关闭
请填处理意见
关闭