24小时热门版块排行榜    

查看: 1133  |  回复: 3
本帖产生 1 个 博学EPI ,点击这里进行查看

闲人与懒猫

木虫 (小有名气)

[求助] matlab 非线性拟合 数据拟合

我在做一个动力学的方程拟合,用的是最小二乘法,下面是程序,不知道哪里有问题,结果偏差很大,求各位大神帮助!
原方程 K=k0*exp(-E/T)*PH^a*PS^b*P,知K T PH PS P,求 k0 E a b
function nihe
global T0 PH0 PS0 P0 K0
T0=[5014.5891 5014.5891 5014.5891 5180.8691 5180.8691 5180.8691 5347.1491 5347.1491 5347.1491];
PH0=[788732.39 592592.59 396039.6 790123.46 394366.19 594059.41 792079.21 591549.29 395061.73];
PS0=[20.19 13.28 7.1 17.71 10.1 10.65 14.2 15.15 8.85];
P0=[24220.91 29411.10 20587.77 20587.77 29411.10 24220.91 29411.10 20587.77 24220.91];
K0=[1047 1122 1233 1188 1225 1148 1815 1185 1125];
y=log(K0);
x1=T0
x2=log(PH0);
x3=log(PS0);
x4=log(P0);
y=y';
X=[ones(size(y)) x1' x2' x3' x4'];
[b bint]=regress(y,X,0.1);
k=exp(b(1));
E=b(2)
a=b(3);
c=b(4);
beta0=[k E a c];
lb=[];ub=[];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
     lsqnonlin(@ObjFunc,beta0,lb,ub);
ci=nlparci(beta,resid,jacobian)

% 残差关于拟合值的残差图
K0c = Rate(beta,PH0,PS0,T0,P0);
plot(K0c,resid,'*')
xlabel('反应速率拟合值, torr s^-^1')
ylabel('残差R, torr s^-^1')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tE = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\ta = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tc = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm)
% ----------------------------------------------------------------
function f=ObjFunc(beta)
global T0 PH0 PS0 P0 K0
f=K0-Rate(beta,T0,PH0,PS0,P0);
% ------------------------------------------------------------------
function K=Rate(beta,T,PH,PS,P)
K=beta(1).*exp(-beta(2)./T).*PH.^beta(3).*PS.^beta(4).*P;



运行结果如下
ci =

  1.0e+004 *

   -0.0000    0.0000
   -2.4270    3.5058
   -0.0000    0.0002
   -0.0002    0.0001

Estimated Parameters:
        k = 0.0001 ± 0.0015
        E = 5393.89 ± 29664.08
        a = 0.64 ± 1.05
        c = 0.64 ± 1.05
        The sum of the squares is: 2.2e+005

» 猜你喜欢

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

最误人二字:功名。
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖置顶 ( 共有1个 )

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
闲人与懒猫: 金币+40, 博学EPI+1, ★★★很有帮助, 谢谢! 2015-01-29 20:14:43
均方差(RMSE):145.616909851711
残差平方和(SSE):190838.559912851
相关系数(R): 0.808436694770782
相关系数之平方(R^2): 0.653569889451906
确定系数(DC): 0.533303922309213
卡方系数(Chi-Square): 84.1949588985453
F统计(F-Statistic): 3.57120350124504

参数        最佳估算
----------        -------------
a        0.000908846493918729
b        15363.6309483451
c        0.613132584682165
d        -0.466403903584844

====== 结果输出 ======

No        实测值y        计算值y
1        1047        1044.65871777826
2        1122        1294.24697276882
3        1233        947.632511671439
4        1188        1042.59182404487
5        1225        1263.95133737708
6        1148        1305.46096405217
7        1815        1813.28448120857
8        1185        1029.71059671882
9        1125        1215.31618898618
2楼2015-01-29 20:07:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

闲人与懒猫

木虫 (小有名气)

引用回帖:
2楼: Originally posted by dingd at 2015-01-29 20:07:38
均方差(RMSE):145.616909851711
残差平方和(SSE):190838.559912851
相关系数(R): 0.808436694770782
相关系数之平方(R^2): 0.653569889451906
确定系数(DC): 0.533303922309213
卡方系数(Chi-Square): 84.1949 ...

谢谢大神,我用1stopt 1.5做了一下,就是这个结果,但是不是想要的数值,忽然发现有个单位是错的,还不知道影响大不大,我再算算。
最误人二字:功名。
3楼2015-01-29 20:16:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

闲人与懒猫

木虫 (小有名气)

单位基本没影响,求助matlab大神。。!
最误人二字:功名。
4楼2015-01-29 22:38:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 闲人与懒猫 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +20 木水思豆 2026-08-25 23/1150 2026-08-26 06:35 by echo8914667
[公派出国] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 PbwVEMK5haaE 2026-08-25 5/250 2026-08-26 05:09 by cNXvBfCpiZOM
[考博] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 Qj98d49dHFb9 2026-08-25 5/250 2026-08-26 04:57 by cNXvBfCpiZOM
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 7K1CJE38xLG4 2026-08-25 4/200 2026-08-26 02:04 by cNXvBfCpiZOM
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了!? +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 19/950 2026-08-25 17:57 by chick875
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
[基金申请] 2026国自然函评费到账 +20 羊腰板 2026-08-21 23/1150 2026-08-25 16:34 by zsna
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +3 苏东坡二世 2026-08-23 3/150 2026-08-25 10:35 by shisan1313
[基金申请] 今天基金会出结果吗?20260819 +17 kkkl_v 2026-08-19 18/900 2026-08-25 09:41 by windflowerwy
[基金申请] 我面上完蛋了 +13 且听虎啸 2026-08-20 14/700 2026-08-25 09:10 by mrkang
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] filecode,4个jtjc了 +14 ziyangfang 2026-08-19 17/850 2026-08-24 18:37 by 哈哈蛤?
[基金申请] 让我中一个面上吧! +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
信息提示
请填处理意见