24小时热门版块排行榜    

查看: 2070  |  回复: 1

wtb1987

新虫 (初入文坛)

[交流] 【求助】MATLAB数值模拟 已有1人参与

我编写了一个程序,结果拟合出来未知参数的置信区间很大,希望大家能帮我指导一下!!!
function HZSM_8
clear all;clc
expdata=[4.28 0.0724 0.0614 0.0075 0.0013 0.0021;
         8.56 0.1126 0.0859 0.0164 0.0039 0.0064;
         16.69 0.1529 0.1034 0.0278 0.0079 0.0138;
         33.76 0.1928 0.1134 0.0412 0.0132 0.0251;
         69.72 0.2305 0.1158 0.0555 0.0192 0.04;];
yexp=expdata(:,2:6);
x0=[0 0 0 0 0];
k0=[1.862 0.578 1.321 0.8 9.813 1.018 0.2409 0.626 4.5301 18.174 40.311 314.1];
options=optimset('Algorithm','Levenberg-Marquardt','LargeScale','off');
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@objfunc,k0,[],[],options,x0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('\t k1=%.4f±%.4f\n',k(1),ci(1,2)-k(1))
fprintf('\t k2=%.4f±%.4f\n',k(2),ci(2,2)-k(2))
fprintf('\t k3=%.4f±%.4f\n',k(3),ci(3,2)-k(3))
fprintf('\t k4=%.4f±%.4f\n',k(4),ci(4,2)-k(4))
fprintf('\t k5=%.4f±%.4f\n',k(5),ci(5,2)-k(5))
fprintf('\t k6=%.4f±%.4f\n',k(6),ci(6,2)-k(6))
fprintf('\t k7=%.4f±%.4f\n',k(7),ci(7,2)-k(7))
fprintf('\t k8=%.4f±%.4f\n',k(8),ci(8,2)-k(8))
fprintf('\t k9=%.4f±%.4f\n',k(9),ci(9,2)-k(9))
fprintf('\t k10=%.4f±%.4f\n',k(10),ci(10,2)-k(10))
fprintf('\t k11=%.4f±%.4f\n',k(11),ci(11,2)-k(11))
fprintf('\t k12=%.4f±%.4f\n',k(12),ci(12,2)-k(12))
fprintf('\t the sum of the squares is:%.1e\n\n',resnorm)
function f=objfunc(k,x0,yexp)
%使用四阶龙格库塔法求取常微分方程组的值
tspan=[0 4.28 8.56 16.69 33.76 69.72];
[t,x]=ode45(@equations,tspan,x0,[],k)
y=x(2:6,;
%目标函数的定义
f1=y(:,1)-yexp(:,1);
f2=y(:,2)-yexp(:,2);
f3=y(:,3)-yexp(:,3);
f4=y(:,4)-yexp(:,4);
f5=y(:,5)-yexp(:,5);
f=[f1;f2;f3;f4;f5];
function dxdt=equations(t,x,k)
P=1;
PT=(1-x(1))*P/(6+x(1));
Ppx=x(2)*P/(6+x(2));
PO=x(4)*P/(6+x(1));
PMx=x(3)*P/(6+x(1));
PTMB=x(5)*P/(6+x(5));
PM=[0.5-(x(1)+x(5))]/(6+x(1));
Z=(1+k(7)*PT+k(8)*PM+k(9)*Ppx+k(10)*PMx+k(11)*PO+k(12)*PTMB)^2;
dx1dt=k(1)*PT*PM/Z;
dx2dt=[k(1)*PT*PM-k(2)*(Ppx-PMx/2.2)-k(4)*Ppx*PM]/Z;
dx3dt=[k(2)*(Ppx-PMx/2.2)-k(3)*(PMx-PO/0.4)-k(5)*PMx*PM]/Z;
dx4dt=[k(3)*(PMx-PO/0.4)-k(6)*PO*PM]/Z;
dx5dt=(k(4)*Ppx*PM+k(5)*PMx*PM+k(6)*PO*PM)/Z;
dxdt=[dx1dt;dx2dt;dx3dt;dx4dt;dx5dt];

结果是:
k1=1.9174±58287.8297
         k2=0.5759±17506.6005
         k3=1.4271±43384.6026
         k4=0.8909±27081.9850
         k5=9.4426±287049.8962
         k6=2.5515±77564.2240
         k7=-0.0804±227820.0664
         k8=1.4310±276004.5168
         k9=5.7924±106309.8775
         k10=17.1832±259330.6709
         k11=37.7388±551649.7103
         k12=314.3976±4536110.8346
         the sum of the squares is:6.2e-007
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

spiderone

木虫 (小有名气)


小木虫: 金币+0.5, 给个红包,谢谢回帖
您这个问题是怎么解决的啊,我的也是这个情况,参数区间大的很!
2楼2013-12-02 15:32:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wtb1987 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[找工作] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 7/350 2026-02-08 07:46 by vs90ilomwc
[公派出国] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 8/400 2026-02-08 07:32 by vs90ilomwc
[考博] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 8/400 2026-02-08 07:27 by vs90ilomwc
[硕博家园] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 8/400 2026-02-08 07:07 by vs90ilomwc
[硕博家园] 博士延得我,科研能力直往上蹿 +8 偏振片 2026-02-02 8/400 2026-02-08 06:52 by liyeqik
[考博] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +3 5lbyq5wrhb 2026-02-07 3/150 2026-02-08 02:52 by vs90ilomwc
[论文投稿] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +3 3rkserf6qr 2026-02-07 4/200 2026-02-08 02:45 by vs90ilomwc
[硕博家园] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +3 3rkserf6qr 2026-02-07 3/150 2026-02-08 02:32 by vs90ilomwc
[硕博家园] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +5 2h7du0nuhk 2026-02-07 5/250 2026-02-08 02:27 by vs90ilomwc
[考博] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +5 2h7du0nuhk 2026-02-07 5/250 2026-02-08 02:25 by vs90ilomwc
[考博] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 6/300 2026-02-08 02:07 by vs90ilomwc
[教师之家] 售SCI一区文章,我:8 O5 51O 54,科目齐全 +4 2h7du0nuhk 2026-02-07 6/300 2026-02-08 02:05 by vs90ilomwc
[教师之家] 有院领导为了换新车,用横向课题经费买了俩车 +7 瞬息宇宙 2026-02-04 7/350 2026-02-07 21:47 by tfang
[有机交流] 酰胺脱乙酰基 10+5 chibby 2026-02-03 12/600 2026-02-07 19:29 by 江东闲人
[基金申请] 同年申请2项不同项目,第1个项目里不写第2个项目的信息,可以吗 +4 hitsdu 2026-02-06 4/200 2026-02-07 13:07 by jurkat.1640
[基金申请] 有时候真觉得大城市人没有县城人甚至个体户幸福 +9 苏东坡二世 2026-02-04 10/500 2026-02-07 12:37 by 小毛球
[考博] 天津大学招2026.09的博士生,欢迎大家推荐交流(博导是本人) +4 a793625982 2026-02-05 5/250 2026-02-07 10:57 by a793625982
[公派出国] CSC & MSCA 博洛尼亚大学能源材料课题组博士/博士后招生|MSCA经费充足、排名优 +4 雨念 2026-02-01 6/300 2026-02-06 23:32 by MelissaPon
[基金申请] 面上项目申报 +3 Tide man 2026-02-01 3/150 2026-02-05 22:56 by god_tian
[教师之家] 遇见不省心的家人很难过 +18 otani 2026-02-03 22/1100 2026-02-04 11:06 by tangmnt
信息提示
请填处理意见