maab数学建模实例.docx
- 文档编号:26103917
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:19
- 大小:18.96KB
maab数学建模实例.docx
《maab数学建模实例.docx》由会员分享,可在线阅读,更多相关《maab数学建模实例.docx(19页珍藏版)》请在冰豆网上搜索。
maab数学建模实例
第四周
3.
functiony=mj()
forx0=0:
:
8
x1=x0^*x0^2+*;
if(abs(x1)<
x0
end
end
4.分别用简单迭代法、埃特金法、牛顿法求解方程,并比较收敛性与收敛速度(分别取10-3、10-5、10-8)。
简单迭代法:
functiony=jddd(x0)
x1=(20+10*x0-2*x0^2-x0^3)/20;
k=1;
while(abs(x1-x0)>=
x0=x1;
x1=(20+10*x0-2*x0^2-x0^3)/20;k=k+1;
end
x1
k
埃特金法:
functiony=etj(x0)
x1=(20-2*x0^2-x0^3)/10;
x2=(20-2*x1^2-x1^3)/10;
x3=x2-(x2-x1)^2/(x2-2*x1+x0);
k=1;
while(abs(x3-x0)>=
x0=x3;
x1=(20-2*x0^2-x0^3)/10;
x2=(20-2*x1^2-x1^3)/10;
x3=x2-(x2-x1)^2/(x2-2*x1+x0);k=k+1;
end
x3
k
牛顿法:
functiony=newton(x0)
x1=x0-fc(x0)/df(x0);
k=1;
while(abs(x1-x0)>=
x0=x1;
x1=x0-fc(x0)/df(x0);k=k+1;
end
x1
k
functiony=fc(x)
y=x^3+2*x^2+10*x-20;
functiony=df(x)
y=3*x^2+4*x+10;
第六周
1.解例6-4(p77)的方程组,分别采用消去法(矩阵分解)、Jacobi迭代法、Seidel迭代法、松弛法求解,并比较收敛速度。
消去法:
x=a\d
或
[L,U]=lu(a);
x=inv(U)inv(L)d
Jacobi迭代法:
functions=jacobi(a,d,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
C=inv(D);
B=C*(L+U);
G=C*d;
s=B*x0+G;
n=1;
whilenorm(s-x0)>=
x0=s;
s=B*x0+G;
n=n+1;
end
n
Seidel迭代法:
functions=seidel(a,d,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
C=inv(D-L);
B=C*U;
G=C*d;
s=B*x0+G;
n=1;
whilenorm(s-x0)>=
x0=s;
s=B*x0+G;
n=n+1;
end
n
松弛法:
functions=loose(a,d,x0,w)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
C=inv(D-w*L);
B=C*((1-w)*D+w*U);
G=w*C*d;
s=B*x0+G;
n=1;
whilenorm(s-x0)>=
x0=s;
s=B*x0+G;
n=n+1;
end
n
2.练习MATLAB的常用矩阵语句,就龙格现象函数(p88)练习插值语句interp,spline,并比较。
3.测得血液中某药物浓度随时间的变化值为:
t(h)
C(mg/L)
求t=,,,时的浓度C.
分别用n=4,5,9的拉格朗日插值计算;并用样条函数插值计算,并比较结果。
拉格朗日插值:
functions=lagr(n)
x=[];
y=[];
x0=[];
m=length(x0);
fori=1:
m
D=abs(x-x0(i));
I=1;
whileI<=n+1
fora=1:
length(x)
ifD(a)==min(D)
c(I)=a;
D(a)=max(D)+1;
break
end
end
I=I+1;
end
b=sort(c);
z=x0(i);
t=;
fork=1:
length(b)
u=;
forj=1:
length(b)
ifj~=k
u=u*(z-x(b(j)))/(x(b(k))-x(b(j)));
end
end
t=t+u*y(b(k));
end
s(i)=t;
end
样条函数差值:
Interp1(x,y,x0,’spline’)
Spline(x,y,x0)
第八周
1.给定某药物浓度随时间的变化值(作业3),1)分别采用样条函数和三点公式(设h=求结点处的导数值,并比较结果。
2)求该时间段的平均浓度(定步长S法)
样条函数:
x=[];
y=[];
pp=csape(x,y,'not-a-knot');
df=fnder(pp);
df1=ppval(df,x)
三点公式:
functiondf=sandian()
t=[];
c=[];
h=;n=length(t);
fori=1:
n
x0=t(i);
y0=c(i);
y1=spline(t,c,x0+h);
y2=spline(t,c,x0+2*h);
y3=spline(t,c,x0-h);
y4=spline(t,c,x0-2*h);
switchi
case1
df(i)=(-3*y0+4*y1-y2)/(2*h);
casen
df(i)=(y4-4*y3+3*y0)/(2*h);
otherwise
df(i)=(y1-y3)/(2*h);
end
end
end
平均浓度:
functionaveragec=simpson()
t=[];
c=[];
m=(t
(1)+t(10))/2;
y=spline(t,c,m);
averagec=(c
(1)+4*y+c(10))/6;
end
2.练习MATLAB常用的trapz,quad,quadl等语句。
计算:
x=0:
8;
y=1./(sqrt(2.*pi)).*exp(-(x-4).^2./2);
z=trapz(x,y)
functiony=jifen(x)
y=1./(sqrt(2.*pi)).*exp(-(x-4).^2./2);
q1=quad('jifen',0,8,
q2=quadl('jifen',0,8,
3.采用变步长经典R-K法,ode23,ode45计算例9-5,并作比较。
变步长经典R-K法:
(可能有问题)
functionz=jdrk(m)
x0=[252]';a=0;b=15;
n=length(x0);
z=zeros(n,m);
k1=zeros(n,1);
k2=zeros(n,1);
k3=zeros(n,1);
k4=zeros(n,1);
t=a;
x=x0;x2=zeros(n,1);x3=x2;x4=x2;
h=choose(m);
m1=15/h+1;
fork=1:
m1
k1=prey(t,x);
fori=1:
n
x2(i)=x(i)+1/2*h*k1(i);
end
k2=prey(t+h/2,x2);
fori=1:
n
x3(i)=x(i)+1/2*h*k2(i);
end
k3=prey(t+h/2,x3);
fori=1:
n
x4(i)=x(i)+h*k3(i);
end
k4=prey(t+h,x4);
fori=1:
n
x(i)=x(i)+h/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i));
z(i,k)=x(i);
end
t=t+h;
end
h1=length(z);
t2=[a:
(b-a)/(h1-1):
b];
plot(t2,z)
gtext('x1(t)')
gtext('x2(t)')
functionh=choose(n)
h=15/(n-1);
t0=0;
x0=[252]';
k11=prey(t0,x0);
k21=prey(t0+h/2,x0+h/2*k11);
k31=prey(t0+h/2,x0+h/2*k21);
k41=prey(t0+h,x0+h*k31);
x1=x0+h/6*(k11+2*k21+2*k31+k41);
k12=prey(t0,x0);
k22=prey(t0+h/4,x0+h/4*k12);
k32=prey(t0+h/4,x0+h/4*k22);
k42=prey(t0+h/2,x0+h/2*k32);
x2=x0+h/12*(k12+2*k22+2*k32+k42);
ifabs(x2-x1)<
whileabs(x2-x1)<
h=h*2;
k11=prey(t0,x0);
k21=prey(t0+h/2,x0+h/2*k11);
k31=prey(t0+h/2,x0+h/2*k21);
k41=prey(t0+h,x0+h*k31);
x1=x0+h/6*(k11+2*k21+2*k31+k41);
k12=prey(t0,x0);
k22=prey(t0+h/4,x0+h/4*k12);
k32=prey(t0+h/4,x0+h/4*k22);
k42=prey(t0+h/2,x0+h/2*k32);
x2=x0+h/12*(k12+2*k22+2*k32+k42);
end
h=h/2;
else
whileabs(x2-x1)>=
h=h/2;
k11=prey(t0,x0);
k21=prey(t0+h/2,x0+h/2*k11);
k31=prey(t0+h/2,x0+h/2*k21);
k41=prey(t0+h,x0+h*k31);
x1=x0+h/6*(k11+2*k21+2*k31+k41);
k12=prey(t0,x0);
k22=prey(t0+h/4,x0+h/4*k12);
k32=prey(t0+h/4,x0+h/4*k22);
k42=prey(t0+h/2,x0+h/2*k32);
x2=x0+h/12*(k12+2*k22+2*k32+k42);
end
end
functionxdot=prey(t,x)
r=1;a=;b=;c=;
xdot=[r-a*x
(2)0;0-b+c*x
(1)]*x;
ode23,ode45:
[t,x]=ode23('prey',[0:
:
15],[252]);
plot(t,x)
[t,x]
gtext('x1(t)')
gtext('x2(t)')
[t,x]=ode45('prey',[0:
:
15],[252]);
plot(t,x)
[t,x]
gtext('x1(t)')
gtext('x2(t)')
第十周
1.熟悉常用的概率分布、概率密度函数图、分位点。
(统计工具箱)
2.对例10-1作统计分组(每组间隔分别为3cm、5cm),并作直方图,计算特征值与置信区间;如假设μ0=175作检验(α=)
functiony=zf(n)
data=[162166171167157168164178170152158153160174159167171168182160159172178166159173161150164175173163165146163162158164169170164179169178170155169160174159168151176164161163172167154164153165161168166166148161163177178171162156165176170156172163165149176170182159164179162151170160165167155168179165184157];
m=ceil((max(data)-min(data))/n);
hist(data,m)
data=[162166171167157168164178170152158153160174159167171168182160159172178166159173161150164175173163165146163162158164169170164179169178170155169160174159168151176164161163172167154164153165161168166166148161163177178171162156165176170156172163165149176170182159164179162151170160165167155168179165184157];
E=mean(data)
D=var(data)
[musigmamucisigmaci]=normfit(data,
[h,p,ci]=ttest(data,175,,0)
3.自行寻找生物学数据,进行分析,试作曲线图、条形图、饼图。
(包括图示)
第十二周
1、作图练习不同形式误差的叠加,随机误差+周期性误差;随机误差+线性误差;随机误差+恒定误差。
(自行设计数据,注意误差数量级的选取)
2、作errorbar图(本课件Page3-A)
T=[];
S=[];
E=[];
errorbar(T,S,E)
xlabel('T/')
ylabel('S/ofwater)')
title('Solubilityof|á-FormGlycineinWater')
3、异常数据剔除
拉依特准则:
functiony=lyt()
x=[];
mu=mean(x);sigma=std(x);n=length(x);
ifn<10
m=2;
elsem=3;
fori=1:
n
ifabs(x(i)-mu)>m*sigma
i
x(i)
end
end
end
格鲁布斯准则:
functiony=grubbs()
x=[];
mu=mean(x);sigma=std(x);n=length(x);
fori=1:
n
ifabs((x(i)-mu)/sigma)>=%格鲁布斯极限值(n=26):
i
x(i)
end
end
end
狄克逊准则:
x=[];
n4=0;
f=[];
whilen4==0
z=sort(x);
n=length(x);
n5=1;
a=(z(3)-z
(1))/(z(n-2)-z
(1));
n1=0;
ifa>f(n5)
n1=1;
z
(1)
end
n2=0;
b=(z(n)-z(n-2))/(z(n)-z(3));
ifb>f(n5)
n2=1;
z(n)
end
x1=[00];
ifn1==1&&n2==0
forn3=1:
(n-1)
x1(n3)=z(n3+1);
end
n5=n5+1;
end
ifn1==1&&n2==1
forn3=1:
(n-2)
x1(n3)=z(n3+1);
end
n5=n5+2;
end
ifn1==0&&n2==1
forn3=1:
(n-1)
x1(n3)=z(n3);
end
n5=n5+1;
end
x=x1;
ifn1==0&&n2==0
n4=1;
end
end
第十四周:
1.大肠杆菌比生长速率测定。
在一定培养条件下,培养大肠杆菌,测得实验数据如下表。
求:
该条件下,大肠杆菌的最大比生长速率μm,半饱和常数Ks,并作模型检验。
S(mg/L)
μ(h-1)
S(mg/L)
μ(h-1)
6
122
13
153
33
170
40
221
64
210
102
s=[613334064102122153170221210];
mu=[];
spmu=s./mu;n=length(s);
a=polyfit(s,spmu,1);
mum=1/a
(1)
ks=a
(2)/a
(1)
lxx=sum(s.^2)-1/n*(sum(s))^2;
lyy=sum(spmu.^2)-1/n*(sum(spmu))^2;
lxy=sum(s.*spmu)-1/n*sum(s)*sum(spmu);
r=lxy/(sqrt(lxx*lyy))
R=corrcoef(s,spmu)
Qr=lxy^2/lxx;
Q=(lxx*lyy-lxy^2)/lxx;
F=Qr/(Q/(n-2))
2.多元线性回归
Pa=[]';
Pb=[]';
Pc=[]';
r=[]';
k0=ones(8,1);alpha=;
r0=log(r);
Pa0=log(Pa);
Pb0=log(Pb);
Pc0=log(Pc);
p=[k0Pa0Pb0Pc0];
[b,bint,r,rint,stats]=regress(r0,p,alpha)
k=exp(b
(1))
m=r'*r
p1=[Pa0Pb0Pc0];
stepwise(p1,r0)
第十六周
1.对作业(7)的两题,分别作非线性回归,并比较参数值和残差。
functiony=ecolinlin(beta0)
s=[613334064102122153170221210]';
mu=[]';
[beta,r,J]=nlinfit(s,mu,'szsl',beta0);
mum=beta
(1)
ks=beta
(2)
r
J
functionmu=szsl(beta,s)
mu=beta
(1).*s./(beta
(2)+s);
end
functiony=reacnlin(beta0)
Pa=[]';
Pb=[]';
Pc=[]';
R=[]';
p=[PaPbPc];
[beta,r,J]=nlinfit(p,R,'fydl',beta0);
k=beta(4)
n1=beta
(1)
n2=beta
(2)
n3=beta(3)
r
J
functionr=fydl(beta,p)
r=beta(4).*p(:
1).^beta
(1).*p(:
2).^beta
(2).*p(:
3).^beta(3);
end
2.作二次正交回归。
数据如右,比较不同模型计算结果。
x1=[1111-1-1-1-10000000000]';
x2=[11-1-111-1-10000000000]';
x3=[1-11-11-11-10000000000]';
y=[]';
x=[x1x2x3];alpha=;
rstool(x,y,'linear',alpha)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- maab 数学 建模 实例