matlab数学建模实例Word文档格式.docx
《matlab数学建模实例Word文档格式.docx》由会员分享，可在线阅读，更多相关《matlab数学建模实例Word文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
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;
n
Seidel迭代法：
functions=seidel（a,d,x0）
C=inv（D-L）;
B=C*U;
=1.0e-5
松弛法：
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
I=I+1;
b=sort（c）;
z=x0（i）;
t=0.0;
fork=1:
length（b）
u=1.0;
forj=1:
ifj~=k
u=u*（z-x（b（j）））/（x（b（k））-x（b（j）））;
t=t+u*y（b（k））;
s（i）=t;
样条函数差值：
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:
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）;
平均浓度：
functionaveragec=simpson（）
m=（t
（1）+t（10））/2;
y=spline（t,c,m）;
averagec=（c
（1）+4*y+c（10））/6;
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）
q1=quad（'
jifen'
0,8,1.0e-8）
q2=quadl（'
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:
x2（i）=x（i）+1/2*h*k1（i）;
k2=prey（t+h/2,x2）;
x3（i）=x（i）+1/2*h*k2（i）;
k3=prey（t+h/2,x3）;
x4（i）=x（i）+h*k3（i）;
k4=prey（t+h,x4）;
x（i）=x（i）+h/6*（k1（i）+2*k2（i）+2*k3（i）+k4（i））;
z（i,k）=x（i）;
t=t+h;
h1=length（z）;
t2=[a:
（b-a）/（h1-1）:
b];
plot（t2,z）
gtext（'
x1（t）'
）
x2（t）'
functionh=choose（n）
h=15/（n-1）;
t0=0;
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）<
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）;
h=h/2;
else
whileabs（x2-x1）>
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]
[t,x]=ode45（'
第十周
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）
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;
ifabs（x（i）-mu）>
m*sigma
i
x（i）
格鲁布斯准则：
functiony=grubbs（）
ifabs（（x（i）-mu）/sigma）>
=2.681%格鲁布斯极限值（n=26）：
0.005-3.1570.01-3.0290.025-2.8410.05-2.681
狄克逊准则：
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）
n2=0;
b=（z（n）-z（n-2））/（z（n）-z（3））;
ifb>
n2=1;
z（n）
x1=[00];
ifn1==1&
&
n2==0
forn3=1:
（n-1）
x1（n3）=z（n3+1）;
n5=n5+1;
n2==1
（n-2）
n5=n5+2;
ifn1==0&
x1（n3）=z（n3）;
x=x1;
n4=1;
第十四周：
1.大肠杆菌比生长速率测定。
在一定培养条件下，培养大肠杆菌，测得实验数据如下表。
求：
该条件下，大肠杆菌的最大比生长速率μm，半饱和常数Ks，并作模型检验。
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）;
functiony=reacnlin（beta0）
R=[1.971.050.730.250.180.130.070.04]'
p=[PaPbPc];
[beta,r,J]=nlinfit（p,R,'
fydl'
k=beta（4）
n1=beta
（1）
n2=beta
（2）
n3=beta（3）
functionr=fydl（beta,p）
r=beta（4）.*p（:
1）.^beta
（1）.*p（:
2）.^beta
（2）.*p（:
3）.^beta（3）;
2.作二次正交回归。
数据如右，比较不同模型计算结果。
x1=[1111-1-1-1-1-1.681.680000000000]'
x2=[11-1-111-1-100-1.681.6800000000]'
x3=[1-11-11-11-10000-1.681.68000000]'
y=[730.2780.5266.7224.5783.1837.5622.6538.3536.2221.2214.2926.2702.4680.1868.5788.3856.5853.4772.6848.4]'
x=[x1x2x3];
rstoo
升级会员 