matlab数学建模实例.docx
- 文档编号:9298898
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:19
- 大小:19.75KB
matlab数学建模实例.docx
《matlab数学建模实例.docx》由会员分享,可在线阅读,更多相关《matlab数学建模实例.docx(19页珍藏版)》请在冰豆网上搜索。
matlab数学建模实例
第四周
3.
functiony=mj()
forx0=0:
0.01:
8
x1=x0^3-11.1*x0^2+38.79*x0-41.769;
if(abs(x1)<1.0e-8)
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)>=1.0e-3)
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)>=1.0e-3)
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)>=1.0e-3)
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)>=1.0e-8
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)>=1.0e-5
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)>=1.0e-8
x0=s;
s=B*x0+G;
n=n+1;
end
n
2.练习MATLAB的常用矩阵语句,就龙格现象函数(p88)练习插值语句interp,spline,并比较。
3.测得血液中某药物浓度随时间的变化值为:
t(h)
0.25
0.5
1.0
1.5
2.0
3.0
4.0
6.0
8.0
10.0
C(mg/L)
19.30
18.15
15.36
14.10
12.89
9.32
7.55
5.24
3.86
2.88
求t=0.45,1.75,5.0,6.0时的浓度C.
分别用n=4,5,9的拉格朗日插值计算;并用样条函数插值计算,并比较结果。
拉格朗日插值:
functions=lagr(n)
x=[0.250.51.01.52.03.04.06.08.010.0];
y=[19.3018.1515.3614.1012.899.327.555.243.862.88];
x0=[0.451.755.06.0];
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=0.0;
fork=1:
length(b)
u=1.0;
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=0.1)求结点处的导数值,并比较结果。
2)求该时间段的平均浓度(定步长S法)
样条函数:
x=[0.250.51.01.52.03.04.06.08.010.0];
y=[19.3018.1515.3614.1012.899.327.555.243.862.88];
pp=csape(x,y,'not-a-knot');
df=fnder(pp);
df1=ppval(df,x)
三点公式:
functiondf=sandian()
t=[0.250.51.01.52.03.04.06.08.010.0];
c=[19.3018.1515.3614.1012.899.327.555.243.862.88];
h=0.1;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=[0.250.51.01.52.03.04.06.08.010.0];
c=[19.3018.1515.3614.1012.899.327.555.243.862.88];
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,1.0e-8)
q2=quadl('jifen',0,8,1.0e-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)<1.0e-5
whileabs(x2-x1)<1.0e-5
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)>=1.0e-5
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=0.1;b=0.5;c=0.02;
xdot=[r-a*x
(2)0;0-b+c*x
(1)]*x;
ode23,ode45:
[t,x]=ode23('prey',[0:
0.1:
15],[252]);
plot(t,x)
[t,x]
gtext('x1(t)')
gtext('x2(t)')
[t,x]=ode45('prey',[0:
0.1:
15],[252]);
plot(t,x)
[t,x]
gtext('x1(t)')
gtext('x2(t)')
第十周
1.熟悉常用的概率分布、概率密度函数图、分位点。
(统计工具箱)
2.对例10-1作统计分组(每组间隔分别为3cm、5cm),并作直方图,计算特征值与置信区间;如假设μ0=175作检验(α=0.05)
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,0.05)
[h,p,ci]=ttest(data,175,0.05,0)
3.自行寻找生物学数据,进行分析,试作曲线图、条形图、饼图。
(包括图示)
第十二周
1、作图练习不同形式误差的叠加,随机误差+周期性误差;随机误差+线性误差;随机误差+恒定误差。
(自行设计数据,注意误差数量级的选取)
2、作errorbar图(本课件Page3-A)
T=[5.012.520.025.028.533.036.046.050.055.0];
S=[141.1166.7198.9226.8241.7259.6283.1334.5354.2384.8];
E=[1.81.50.71.50.20.51.21.11.21.5];
errorbar(T,S,E)
xlabel('T/?
?
')
ylabel('S/(g.kg-1ofwater)')
title('Solubilityof|á-FormGlycineinWater')
3、异常数据剔除
拉依特准则:
functiony=lyt()
x=[25.30725.11225.32425.30025.29525.29325.29425.31425.34125.31525.31425.29925.30325.31325.31125.59025.30925.31625.31025.31725.30625.29125.32525.01025.31525.438];
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=[25.30725.11225.32425.30025.29525.29325.29425.31425.34125.31525.31425.29925.30325.31325.31125.59025.30925.31625.31025.31725.30625.29125.32525.01025.31525.438];
mu=mean(x);sigma=std(x);n=length(x);
fori=1:
n
ifabs((x(i)-mu)/sigma)>=2.681%格鲁布斯极限值(n=26):
0.005-3.1570.01-3.0290.025-2.8410.05-2.681
i
x(i)
end
end
end
狄克逊准则:
x=[25.30725.11225.32425.30025.29525.29325.29425.31425.34125.31525.31425.29925.30325.31325.31125.59025.30925.31625.31025.31725.30625.29125.32525.01025.31525.438];
n4=0;
f=[0.3990.4060.4130.4210.4300.4400.4500.4620.4750.4900.5070.5250.546];
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
0.06
122
0.60
13
0.12
153
0.66
33
0.24
170
0.69
40
0.31
221
0.70
64
0.43
210
0.73
102
0.53
s=[613334064102122153170221210];
mu=[0.060.120.240.310.430.530.600.660.690.700.73];
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=[9.08.68.47.57.06.86.56.0]';
Pb=[8.37.06.24.23.93.52.62.2]';
Pc=[2.74.45.48.39.19.710.911.8]';
r=[1.971.050.730.250.180.130.070.04]';
k0=ones(8,1);alpha=0.05;
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=[0.060.120.240.310.430.530.600.660.690.700.73]';
[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=[9.08.68.47.57.06.86.56.0]';
Pb=[8.37.06.24.23.93.52.62.2]';
Pc=[2.74.45.48.39.1
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 数学 建模 实例
![提示](https://static.bdocx.com/images/bang_tan.gif)