24小时热门版块排行榜    

查看: 2681  |  回复: 6

lxyy

银虫 (小有名气)

[求助] matlab拟合拟均相动力学参数

模型如附件所示,要求里面的k+,k-,m,n,p,我自己编了个程序,但是一直有问题
M文件
function f=fortyfour_D(c,t,rD)
t=[5,10,15,20,30,45,60,80,100,120,150,180,210,240,300,360];
CD=[749.16 ,994.49 ,1382.00,1589.18,1909.81,2222.86,2525.68,2721.90,2896.67,2951.42,2960.67,2999.60,3017.97,3048.91,3072.27,3063.85];
knots = 3;K = 3;      %三次B样条
sp=spap2(knots,K,t,CD);
sp=spap2(newknt(sp),K,t,CD);
pp=fnder(sp);          % 计算B样条函数的导函数
dCDdt=fnval(pp,t);     % 计算t处的导函数值
Wcat=85.6312   %催化剂浓度
rD=dCDdt./Wcat;
f=c(1)*CA.^c(2)*CB.^c(3)-c(4)*CD.^c(5)))-rD;
命令如下
c0=[1,1,1,1,1];
for i=1:100
c=lsqnonlin('fortyfour_D ',c0);
c0=c;
end
c
% 绘制图形
ti = linspace(t(1),t(end),200);    %y=linspace(a,b,n)生成一个行向量,该向量将a与b之间平分为n个点,包含端点a和b。
CDi = fnval(sp,ti); %计算样条函数区间内任意一点的值
plot(t,CD,'ro',ti,CDi,'b-'),xlabel('t'),ylabel('C_D')
legend('实验数据','拟合曲线')
c=[1,1,1,1,1];for i=1:100;c=lsqnonlin('fortyfour_D ',c);c;end
CD=[749.16 ,994.49 ,1382.00,1589.18,1909.81,2222.86,2525.68,2721.90,2896.67,2951.42,2960.67,2999.60,3017.97,3048.91,3072.27,3063.85];
CA = [2732.43 2468.85 2223.97 2023.64 1696.41 1343.42 1054.34 902.20 756.91 717.86 703.53 631.00 652.41 649.53 624.92 744.95];      
CB= [3127.04 2722.48 2414.51 2267.55 1927.70 1500.17 1285.07 1111.91 924.27 886.32 830.86 821.82 801.97 780.39 794.41 783.73];      
rD=[0.5631,0.5282,0.4934,0.4585,0.3888,0.2842,0.1796,0.0482,0.0404,0.0326,0.0209,0.0092,0.0046,0.0040,0.0028,0.0015];
>> c0=[1,1,1,1,1];
for i=1:1000
c=lsqcurvefit('ffD ',c0,CD,rD);
c0=c;
end
c

请大家帮忙看看我程序哪里有问题[ Last edited by lxyy on 2011-12-30 at 16:05 ]
回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : clip_image002.gif
  • 2011-12-30 16:05:36, 503 bytes

» 猜你喜欢

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

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

dingd

铁杆木虫 (职业作家)


dbb627(金币+1): 欢迎交流 2011-12-31 08:41:48
不知理解的对否,用1stOpt试下:
CODE:
Variable CD,CA,CB,rD;
Function rD=k1*Ca^m*Cb^n-k2*Cd^p;
Data;
CD=[749.16,994.49,1382.00,1589.18,1909.81,2222.86,2525.68,2721.90,2896.67,2951.42,2960.67,2999.60,3017.97,3048.91,3072.27,3063.85];
CA=[2732.43,2468.85,2223.97,2023.64,1696.41,1343.42,1054.34,902.20,756.91,717.86,703.53,631.00,652.41,649.53,624.92,744.95];
CB=[3127.04,2722.48,2414.51,2267.55,1927.70,1500.17,1285.07,1111.91,924.27,886.32,830.86,821.82,801.97,780.39,794.41,783.73];
rD=[0.5631,0.5282,0.4934,0.4585,0.3888,0.2842,0.1796,0.0482,0.0404,0.0326,0.0209,0.0092,0.0046,0.0040,0.0028,0.0015];

均方差(RMSE): 0.0187338627500751
残差平方和(SSE): 0.00561532181661841
相关系数(R): 0.996169068786095
相关系数之平方(R^2): 0.992352813606156
决定系数(DC): 0.99235277918993

参数        最佳估算
----------        -------------
k1        14.3465921423275
m        0.0595653946350665
n        -0.452413851540345
k2        1.79201130280021E-8
p        2.22793339858939


2楼2011-12-30 22:02:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
: Originally posted by dingd at 2011-12-30 22:02:01:
不知理解的对否,用1stOpt试下:
[code]
Variable CD,CA,CB,rD;
Function rD=k1*Ca^m*Cb^n-k2*Cd^p;
Data;
CD=[749.16,994.49,1382.00,1589.18,1909.81,2222.86,2525.68,2721.90,2896.67,2951.42,2960.67,29 ...

作为一个菜鸟来说我不知道你这样算对不对,但是以我的经验来看mn应该接近1,p算出来应该接近2才对
3楼2011-12-30 22:39:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
lxyy(金币+1): 有帮助 2012-01-05 16:34:35
如果公式和数据没错,参数也没有范围限制的话,2楼给出的结果从数学角度应该是最好的了,可以将该结果和你知道的结果代入验证一下就知道了。
4楼2011-12-31 08:18:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

M文件
function f=q(c,t,CD,rD)
t=[0 0.08333 0.166667 0.25 0.5 0.75 1 1.5 2 2.666667 3.333333 4 5];
CD=[0 0.53259 1.154323 1.621033 2.644088 3.163488 3.461446 3.710468 3.784974 3.746928 3.761759 3.713404 3.643825];
rD=[0.0710291014584411 0.0700187830917774 0.0655437729479848 0.0576582625377526 0.0329826401194250 0.0166083513736405 0.00912074378901541 0.00278932296049053 0.000202497370868135 0 0 0 0];
f= c(1)*((4.6026-CD).^c(2)*(5.064335-CD).^c(3)-c(4)*CD.^c(5))-rD;
命令如下
c0=[0.1,0.1,0.1,0.1,0.1];
for i=1:50
c=lsqnonlin('q ',c0);
c0=c;
end
c
我把程序简化了下,还是一直运行不了
5楼2011-12-31 11:56:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bao2007pin

铜虫 (正式写手)

表示我不懂!坐等高人吧!!!!!!
6楼2011-12-31 12:43:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxyy

银虫 (小有名气)

引用回帖:
: Originally posted by bao2007pin at 2011-12-31 12:43:36:
表示我不懂!坐等高人吧!!!!!!

谢谢
7楼2011-12-31 14:00:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 lxyy 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 基金啊基金 +3 longfie172 2026-08-20 3/150 2026-08-21 00:04 by wuchongjun
[基金申请] 科研孤儿太难了 +14 我4大白菜 2026-08-20 15/750 2026-08-20 22:33 by 袁向阳007
[基金申请] 让我中一个面上吧! +8 大萍1987 2026-08-20 10/500 2026-08-20 21:48 by 大萍1987
[基金申请] 言之凿凿,明天官网会通知放榜? +15 医学老男孩 2026-08-20 17/850 2026-08-20 20:35 by CaoLG
[基金申请] 时间戳又变了 +12 wuchongjun 2026-08-20 17/850 2026-08-20 20:05 by spzhao
[基金申请] 今天系统多次维护,明天很可能放榜! +10 zju2000 2026-08-16 11/550 2026-08-20 20:02 by cl479861084
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] 只有每年这种时候来逛逛小木虫 +21 yaoyewhu2008 2026-08-20 22/1100 2026-08-20 17:03 by yaoyewhu2008
[基金申请] 放榜前的不淡定 40+3 snowwithsea 2026-08-19 12/600 2026-08-20 15:56 by lfy8008
[基金申请] 重要消息,中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[基金申请] filecode,4个jtjc了 +11 ziyangfang 2026-08-19 13/650 2026-08-19 21:28 by aasahr
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 快农历七夕节了,轻松一下,男人悄悄话,女施主请不要进来。 +6 Tide man 2026-08-14 7/350 2026-08-19 09:56 by ZJTJZ
[基金申请] 什么时候开奖? +6 CrisMessi 2026-08-18 6/300 2026-08-19 08:06 by Equinoxhua
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 各位道友,我要去昆明玩几天,回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
信息提示
请填处理意见