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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab拟合方程参数时初值的选择
南方科技大学公共卫生及应急管理学院2025级博士研究生招生报考通知
161/2返回列表12下一页
查看: 3401  |  回复: 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下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭