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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 动力学方程组参数求解问题
101/1返回列表
查看: 2578  |  回复: 9
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

最自私的女人

新虫 (初入文坛)

[求助] 动力学方程组参数求解问题 已有3人参与

本人动力学方程如下：
dx1/dt=k21x2+k31x3+k41x4+k51x5+k61x6+k71x7+k81x8
dx2/dt=-k21x2+k32x3+k42x4+k52x5+k62x6+k72x7+k82x8
dx3/dt=-k31x3-k32x3+k43x4+ k53x5+k63x6+k73x7+k83x8
dx4/dt=-k41x4-k42x4-k43x4+ k54x5+k64x6+k74x7+k84x8
dx5/dt=-k51x5-k52x5-k53x5-k54x5-k56x52-k57x52+k65x6+k75x7+k85x8
dx6/dt=-k61x6-k62x6-k63x6-k64x6-k65x6-k67x62+k56x52+k76x7+k86x8
dx7/dt=-k71x7-k72x7-k73x7-k74x7-k75x7-k76x7-k78x72+k57x52+k67x62+k87x8
dx8/dt=-k81x8-k82x8-k83x8-k84x8-k85x8-k86x8-k87x8-k89x82+k78x72
dx9/dt=k89x82
方程中有多个缩合反应为表观二级（即平方），已上传图片供参考
k为动力学参数，k大于0，共33个
x为质量分数，x＞0，九个x之和为1
已知数据如下：
温度        时间                    x1            x2           x3                   x4                     x5               x6             x7            x8            x9
430        0.00         0.0192         0.0236         0.1006         0.3284         0.1647         0.1327         0.1981         0.0319         0.0008
430        20.00         0.0287         0.0503         0.1313         0.3156         0.1391         0.1206         0.1755         0.0366         0.0022
430        40.00         0.0349         0.0858         0.1827         0.3058         0.1203         0.0805         0.1476         0.0384         0.0039
430        60.00         0.0423         0.1250         0.2086         0.2950         0.1005         0.0624         0.1180         0.0428         0.0053
本人想编写一个MATLAB程序实现以下功能：
设k初值均为0，求解方程中的33个参数，将求得的参数代入方程中，根据t=0min的数据求得t=20、40、60min的x计算值，将x计算值与x实验值比较，如果误差＜5%，得到的参数就是最终参数，如果误差大于5%，就将得到的参数最为初值        循环上述计算，直到        误差＜5%。       
哪位大侠能帮忙编写一下，小女子不胜感激呢~

动力学方程组参数求解问题
动力学方程组.jpg
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-09-02 13:55:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
参数估计问题或者叫常微分方程组参数拟合。参数太多，MATLAB求解估计困难不小，而且数据仅有4组。
建议你找有高版本1stOpt软件的专家帮你做一下，推荐本版块的专家dingd。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
2楼2014-09-02 14:26:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Cle_key

铁虫 (小有名气)

【答案】应助回帖

感谢参与，应助指数 +1
虽然楼主自称小女子。。。但肯定不帮

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼
坚持就能成功
3楼2014-09-02 17:38:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

polypro

木虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
按自由约束条件，具有唯一解的情况是 自由变量数=组分数-平衡相数+2，LZ的当前数据应该有很多不定解
赞一下 回复此楼
泉涸,鱼相与处于陆,相呴以湿,相濡以沫,不如相忘于江湖。
4楼2014-09-03 08:20:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

最自私的女人

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 月只蓝 at 2014-09-02 14:26:53
参数估计问题或者叫常微分方程组参数拟合。参数太多，MATLAB求解估计困难不小，而且数据仅有4组。
建议你找有高版本1stOpt软件的专家帮你做一下，推荐本版块的专家dingd。

因为论文需要得用MATLAB的，我也知道很难，有很多解，只是想得到个稍微好点的解，我就是不知道判断误差那一部分该怎么编写，我现在一个能稍微计算下的程序呢
赞一下 回复此楼
5楼2014-09-03 10:16:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

最自私的女人

新虫 (初入文坛)
引用回帖:
3楼: Originally posted by Cle_key at 2014-09-02 17:38:33
虽然楼主自称小女子。。。但肯定不帮

赠人玫瑰，手有余香的亲
赞一下 回复此楼
6楼2014-09-03 10:17:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

最自私的女人

新虫 (初入文坛)
引用回帖:
4楼: Originally posted by polypro at 2014-09-03 08:20:26
按自由约束条件，具有唯一解的情况是 自由变量数=组分数-平衡相数+2，LZ的当前数据应该有很多不定解

我也知道有很多不定解呢，我想在这些不定解中找到使得x计算值与x实验值相差最小的那一组，请问你会不会编写误差判断这一部分啊？我自己有个能计算的程序呢
赞一下 回复此楼
7楼2014-09-03 10:18:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

最自私的女人

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 月只蓝 at 2014-09-02 14:26:53
参数估计问题或者叫常微分方程组参数拟合。参数太多，MATLAB求解估计困难不小，而且数据仅有4组。
建议你找有高版本1stOpt软件的专家帮你做一下，推荐本版块的专家dingd。

你好，能麻烦你帮我看一下我的下面这个程序为什么运行不了嘛？
function fourLumpK_VR440
clc;
clear all;
format long;
fourLumpData=[...
0        69.33        21.33        7.66        1.76
15        52.96        26.45        17.84        3.18
30        49.9        27.61        18.29        3.92
45        39.03        28.51        27.93        5
60        33.45        25.4        34.49        7.1
90        25.76        24.99        42.38        7.48
];
% vector of yields/%
Xexp=fourLumpData(2:6,2:5);
% oil residence time,unit: min
time= fourLumpData(2:6,1)';  

% initiate parameter for optimization
X0=fourLumpData(1,2:5);
%      k1      k2     k3      k4     k5  
K0 =[0.0014 0.0012 0.0005 0.0020 0.0013];
lb=0;
ub=inf;

%计时开始拟合
tic;

K=lsqnonlin(@objFun,K0,lb,ub,[],X0,Xexp,time);
% 拟合结果显示
disp(K)
% 计时结束
disp('耗时：');
toc;
% 保存拟合结果
save fourLumpK_VR440;

% 检验-----------------------------------------------
tf=max(time);
tspan=[0:1:tf];
[nr,nc]=size(Xexp);
XcalFinal=zeros(nr,nc);
for i=1:nr
        [t,Xcal]=ode45(@modelEquation,tspan,X0,[],K,tf);
    n=find(t>=time(i));
        m=n(1);
        %线性内插
        Xc=Xcal(m,-(Xcal(m,-Xcal(m-1,)/(t(m)-t(m-1))*(t(m)-time(i));
        XcalFinal(i,=Xc;
    save XcalFinal
end
objX = abs((XcalFinal - Xexp)./100).^2;
F=sum(sum(objX))^0.5;
error=(XcalFinal - Xexp)
avrError=sum(abs(error))/nr
relatError=(XcalFinal - Xexp)./Xexp*100
avrRelError=sum(abs(relatError))/nr

rou_2=1-sum((XcalFinal - Xexp).^2)./sum(Xexp.^2)
Fc_11=(sum(XcalFinal.^2)-sum((XcalFinal - Xexp).^2))/5./(sum((XcalFinal - Xexp).^2)/(5*4-5))

% 残差图-----------------------------------------------------------
figure(1);
clf;
hold on;
plot(XcalFinal(:,1),Xexp(:,1),'m*');
plot(XcalFinal(:,2),Xexp(:,2),'r*');
plot(XcalFinal(:,3),Xexp(:,3),'b*');
plot(XcalFinal(:,4),Xexp(:,4),'k*');
legend('尾油','蜡油','轻油','气+焦'); % 图形注解
xlabel('Xcal,%'); % x轴注解
ylabel('Xexp,%'); % y轴注解
% title('Xexp-Xcal'); % 图形标题
hold on;
x=0:0.1:80;
plot(x,x,'k-');

% y-t图形显示
figure(2);
clf;
hold on;
tf=max(time);
tspan=[0:1:tf];
[t,XcalFinal]=ode45(@modelEquation,tspan,X0,[],K,tf);
plot(t,XcalFinal(:,1),'m');
plot(t,XcalFinal(:,2),'r');
plot(t,XcalFinal(:,3),'b');
plot(t,XcalFinal(:,4),'k');
legend('尾油','蜡油','轻油','气+焦'); % 图形注解
xlabel('t/min'); % x轴注解
ylabel('X/%'); % y轴注解
% title('X--t'); % 图形标题
hold;

%绘制试验点
figure(2);
hold on;
plot(fourLumpData(:,1)',fourLumpData(:,2),'m*');
plot(fourLumpData(:,1)',fourLumpData(:,3),'r*');
plot(fourLumpData(:,1)',fourLumpData(:,4),'b*');
plot(fourLumpData(:,1)',fourLumpData(:,5),'k*');

% -----------------------------------------------------------------
function F = objFun(K,X0,Xexp,time)
tf=max(time);
tspan=[0:1:tf];
[nr,nc]=size(Xexp);
XcalFinal=zeros(nr,nc);
for i=1:nr
        [t,Xcal]=ode45(@modelEquation,tspan,X0,[],K,tf);
    n=find(t>=time(i));
        m=n(1);
        %线性内插
        Xc=Xcal(m,-(Xcal(m,-Xcal(m-1,)/(t(m)-t(m-1))*(t(m)-time(i));
        XcalFinal(i,=Xc;
end
objX = abs((XcalFinal - Xexp)./100).^2;
F=sum(sum(objX))^0.5;
error=(XcalFinal - Xexp)
avrError=sum(abs(error))/nr
avrErr=sum(avrError)/nc
relatError=(XcalFinal - Xexp)./Xexp*100
avrRelError=sum(abs(relatError))/nr
avrRelErr=sum(avrRelError)/nc
K

function dxdt = modelEquation(t,X,K,tf)
f1=-(K(1)+K(2)+K(3))*X(1);
f2= K(1)*X(1)-(K(4)+K(5))*X(2);
f3= K(2)*X(1)+K(4)*X(2);
f4= K(3)*X(1)+K(5)*X(2);
dxdt = [f1;f2;f3;f4];
赞一下 回复此楼
8楼2014-09-04 20:20:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

最自私的女人

新虫 (初入文坛)
引用回帖:
8楼: Originally posted by 最自私的女人 at 2014-09-04 20:20:33
你好，能麻烦你帮我看一下我的下面这个程序为什么运行不了嘛？
function fourLumpK_VR440
clc;
clear all;
format long;
fourLumpData=;
% vector of yields/%
Xexp=fourLumpData(2:6,2:5);
% oil residen ...

发现有乱码哎，我上传个附件给你，麻烦你帮我看一下吧。
赞一下 回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : fourLumpK_VR440.m
    • 2014-09-04 20:22:37, 3.21 K
9楼2014-09-04 20:22:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

恋雨zz

新虫 (初入文坛)
有人解出了么
回复此楼
10楼2015-05-02 15:10:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 最自私的女人 的主题更新
101/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[公派出国] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 k0dTPqJtl0jt 2026-08-14 4/200 2026-08-15 04:45 by 4wMiSEwB6436
[考博] 售SCI文章，我:8O.5.5.1O.54,科目全,可十急 +4 k0dTPqJtl0jt 2026-08-14 5/250 2026-08-15 04:45 by 4wMiSEwB6436
[公派出国] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +3 7lpolszZVXgi 2026-08-14 5/250 2026-08-15 04:41 by 4wMiSEwB6436
[找工作] 售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急 +4 k0dTPqJtl0jt 2026-08-14 4/200 2026-08-15 02:21 by 4wMiSEwB6436
[基金申请] 小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊 +11 Tide man 2026-08-10 12/600 2026-08-15 02:12 by home3163
[基金申请] 各位道友，我要去昆明玩几天，回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
[考研] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +3 HFw0lei2R37i 2026-08-14 7/350 2026-08-15 00:40 by 4wMiSEwB6436
[论文投稿] 售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急 +3 HFw0lei2R37i 2026-08-14 5/250 2026-08-15 00:16 by 4wMiSEwB6436
[教师之家] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 HFw0lei2R37i 2026-08-14 5/250 2026-08-14 22:32 by 4wMiSEwB6436
[基金申请] 欢迎发来filecode的Mz6后的代码验证其规律 +23 医学老男孩 2026-08-13 49/2450 2026-08-14 20:02 by zhaifei
[基金申请] 哪位老哥知道今年的国自然具体哪一天放榜？ +5 Ldrop2023 2026-08-13 5/250 2026-08-14 18:48 by ssxclkj
[基金申请] 奇怪，两个人的filecode固定段从头到尾一模一样 +8 布布和一二 2026-08-10 11/550 2026-08-14 14:58 by Equinoxhua
[基金申请] 有时候，自然基金真的不能太认真 （我的申报经验） +6 majunge000 2026-08-11 8/400 2026-08-14 14:44 by 医学老男孩
[硕博家园] 请教兼职经验 +3 是阿文鸭 2026-08-09 3/150 2026-08-14 12:07 by HER12025
[硕博家园] 读博的好处 +4 lnee 2026-08-11 4/200 2026-08-14 10:20 by ahsoarli
[文学芳草园] 阿姨 +4 汪汪锅 2026-08-09 4/200 2026-08-13 19:43 by arzu_hma
[基金申请] 不应该看fileCode +7 且听虎啸 2026-08-12 9/450 2026-08-13 14:27 by flydreamws
[基金申请] 结合人工智能，周易传统文化，filecode打分制来了，3分以上希望很大。 +3 Tide man 2026-08-12 4/200 2026-08-13 08:35 by ZJTJZ
[基金申请] 2019年青年基金涵评意见，大家看看几个A，几个B? +11 Tide man 2026-08-11 11/550 2026-08-13 07:35 by 撸猫猫
[基金申请] 为什么网上很多人说本周 12号出结果 +6 瞬息宇宙 2026-08-10 7/350 2026-08-11 19:25 by Tide man
信息提示
关闭
请填处理意见
关闭