matlab能做数值模拟,【求助】MATLAB数值模拟

我编写了一个程序,结果拟合出来未知参数的置信区间很大,希望大家能帮我指导一下!!!

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,

smile.gif;

%目标函数的定义

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值