matlab数学建模实例0001

1、3.fun cti on y=mj()for x0=0:0.01:8x1=x0A3-11.1*x0A2+38.79*x0-41.769;if (abs(x1)=1.0e-3)x0=x1;x1= (20+10*x0-2*x0A2-x0A3)/20;k=k+1;endx1 k埃特金法：fun cti ony=etj(x0)x1=(20-2*x0A2-x0A3)/10;x2=(20-2*x1A2-x1A3)/10;x3=x2-(x2-x1)A2/(x2-2*x1+x0);k=1;while (abs(x3-x0)=1.0e-3)x0=x3;x1=(20-2*x0A2-x0A3)/10;x2=(20-

2、2*x1A2-x1A3)/10;x3=x2-(x2-x1)A2/(x2-2*x1+x0);k=k+1;endx3k牛顿法：fun cti ony=n ewt on( 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;x1kfunctiony=fc(x)y=xA3+2*xA2+10*x-20;functiony=df(x)y=3*xA2+4*x+10;第六周1. 解例6-4 (p77) 的方程组，分别采用消去法(矩阵分解) 、Jacobi 迭代法、 Seidel 迭代法、松弛 法

3、求解，并比拟收敛速度。消去法：x=ad或L,U=lu(a);x=inv(U)inv(L)dJacobi 迭代法：function s=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;while norm(s-x0)=1.0e-8x0=s;s=B*x0+G;n=n+1;endnSeidel 迭代法：function s=seidel(a,d,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);C=inv(D-L);B=C

4、*U;G=C*d;s=B*x0+G;n=1;while norm(s-x0)=1.0e-5x0=s;s=B*x0+G;n=n+1;end松弛法：fun cti on s=loose(a,d,xO,w) D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);C=i nv(D-w*L);B=C*(1-w)*D+w*U);G=w*C*d;s=B*xO+G;n=1;whilen orm(s-x0)=1.0e-8x0=s;s=B*xO+G;n=n+1;endn2. 练习MATLAB的常用矩阵语句，就龙格现象函数p88 练习插值语句in terp,spl ine ,并比拟。3

5、. 测得血液中某药物浓度随时间的变化值为：t(h)0.250.51.01.52.03.04.06.08.010.0C(mg /L)19.3018.1515.3614.1012.899.327.555.243.862.88求t=0.45, 1.75, 5.0, 6.0时的浓度 C.分别用n=4,5,9的拉格朗日插值计算；并用样条函数插值计算，并比拟结果拉格朗日插值:fun cti ons=lagr (n)x=0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0;y=19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2

6、.88;x0=0.45 1.75 5.0 6.0;m=le ngth(x0);for i=1:mD=abs(x-x0(i);I=1;while l=n+1for a=1:le ngth(x)if D( a)=mi n(D)c(I)=a;D(a)=max(D)+1;breakendI=I+1;endb=sort(c);z=x0(i);t=0.0;for k=1:length(b)u=1.0;for j=1:length(b)if j=ku=u*(z-x(b(j)/(x(b(k)-x(b(j);endendt=t+u*y(b(k);ends(i)=t;end样条函数差值：Interp1(x,y,x

7、0,spline )Spline(x,y,x0)第八周h=0.1)1. 给定某药物浓度随时间的变化值 ( 作业 3) ，1)分别采用样条函数和三点公式求结点处的导数值，并比拟结果。2)求该时间段的平均浓度(定步长S法)样条函数：x=0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0;y=19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88;pp=csape(x,y, not-a-knot );df=fnder(pp);df1=ppval(df,x)三点公式：function df=sandian()t=0.

8、25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0;c=19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88;h=0.1;n=length(t);for i=1:nx0=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);switch icase 1df(i)=(-3*y0+4*y1-y2)/(2*h);case ndf(i)=(y4-4*y3+3*y0)/(2*h);

9、otherwisedf(i)=(y1-y3)/(2*h);endendend平均浓度：function averagec=simpson()t=0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0;c=19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88; m=(t(1)+t(10)/2;y=spline(t,c,m); averagec=(c(1)+4*y+c(10)/6;end2. 练习MATLAB常用的trapz, quad, quadl等语句。计算：x=0:8;y=1./(sqrt(2.*pi).*e

10、xp(-(x-4).A2./2);z=trapz(x,y)function y=jifen(x)y=1./(sqrt(2.*pi).*exp(-(x-4).A2./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法：(可能有问题)function z=jdrk(m)x0=25 2;a=0;b=15;n=length(x0);z=zeros(n,m);k1=zeros(n,1);k2=zeros(n,1);k3=zeros(n,1);k

11、4=zeros(n,1);t=a;x=x0;x2=zeros(n,1);x3=x2;x4=x2;h=choose(m);m1=15/h+1;for k=1:m1k1=prey(t,x);for i=1:nx2(i)=x(i)+1/2*h*k1(i);endk2=prey(t+h/2,x2);for i=1:nx3(i)=x(i)+1/2*h*k2(i);end k3=prey(t+h/2,x3);for i=1:n x4(i)=x(i)+h*k3(i); endk4=prey(t+h,x4);for i=1:n x(i)=x(i)+h/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i

12、); z(i,k)=x(i);endt=t+h;end h1=length(z);t2=a:(b-a)/(h1-1):b; plot(t2,z) gtext( x1(t)gtext( x2(t)function h=choose(n) h=15/(n-1);t0=0; x0=25 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=pr

13、ey(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);if abs(x2-x1)1.0e-5while abs(x2-x1)=1.0e-5h=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(t

14、0,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); endendfunction xdot=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,0:0.1:15,25 2);,0:0.1:15,25 2);t,x=ode23( prey plot(t,x)t,xgtext( x1(t)gtext( x2

15、(t)t,x=ode45( prey plot(t,x)t,xgtext( x1(t)gtext( x2(t)第十周1熟悉常用的概率分布、概率密度函数图、分位点。统计工具箱2对例 10-1 作统计分组 ( 每组间隔分别为 3cm、 5cm) ，并作直方图，计算特征值与置信区 间；如假设 g=175作检验(0=0.05 )function y=zf(n)data=162 166 171 167 157 168 164 178 170 152 158 153 160 174 159 167 171 168 182160 159 172 178 166 159 173 161 150 164 175

16、 173 163 165 146 163 162 158 164 169 170164 179 169 178 170 155 169 160 174 159 168 151 176 164 161 163 172 167 154 164 153165 161 168 166 166 148 161 163 177 178 171 162 156 165 176 170 156 172 163 165 149176 170 182 159 164 179 162 151 170 160 165 167 155 168 179 165 184 157; m=ceil(max(data)-min(

17、data)/n);hist(data,m)data=162 166 171 167 157 168 164 178 170 152 158 153 160 174 159 167 171 168 182160 159 172 178 166 159 173 161 150 164 175 173 163 165 146 163 162 158 164 169 170164 179 169 178 170 155 169 160 174 159 168 151 176 164 161 163 172 167 154 164 153165 161 168 166 166 148 161 163 1

18、77 178 171 162 156 165 176 170 156 172 163 165 149 176 170 182 159 164 179 162 151 170 160 165 167 155 168 179 165 184 157; E=mean(data)D=var(data)mu sigma muci sigmaci=normfit(data,0.05)h,p,ci=ttest(data,175,0.05,0)3自行寻找生物学数据，进行分析，试作曲线图、条形图、饼图。包括图示 第十二周1、作图练习不同形式误差的叠加，随机误差 + 周期性误差；随机误差 + 线性误差；随机误差恒

19、定误差。 自行设计数据，注意误差数量级的选取2、作 errorbar 图(本课件 Page 3-A )T=5.0 12.5 20.0 25.0 28.5 33.0 36.0 46.0 50.0 55.0;S=141.1 166.7 198.9 226.8 241.7 259.6 283.1 334.5 354.2 384.8;E=1.8 1.5 0.7 1.5 0.2 0.5 1.2 1.1 1.2 1.5;errorbar(T,S,E)xlabel( T/? )ylabel( S/(g.kg-1 of water) )title(Solubility of |a -Form Glycine

20、in Water)3、异常数据剔除拉依特准那么：function y=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.315 25.438;mu=mean(x);sigma=std(x);n=length(x);if nm*sigmaix(i) end end end格鲁布斯准那么：function y=grubbs()25.31525.

21、31425.299x=25.30725.112 25.324 25.300 25.295 25.293 25.294 25.314 25.34125.303 25.313 25.311 25.590 25.309 25.316 25.310 25.317 25.306 25.291 25.325 25.01025.315 25.438; mu=mean(x);sigma=std(x);n=length(x);for i=1:nif abs(x(i)-mu)/sigma)=2.681%格鲁布斯极限值( n=26 )： i x(i) endendend狄克逊准那么：25.29925.010x=25

22、.307 25.112 25.324 25.300 25.295 25.293 25.294 25.314 25.341 25.315 25.31425.303 25.313 25.311 25.590 25.309 25.316 25.310 25.317 25.306 25.291 25.32525.315 25.438;n4=0;f=0.399 0.406 0.413 0.421 0.430 0.440 0.450 0.462 0.475 0.490 0.507 0.525 0.546;while n4=0z=sort(x);n=length(x);n5=1;a=(z(3)-z(1)/(

23、z(n-2)-z(1);n1=0;if af(n5)n1=1;z(1)endn2=0; b=(z(n)-z(n-2)/(z(n)-z(3);if bf(n5)n2=1;z(n) endx1=0 0;if n1=1&n2=0for n3=1:(n-1) x1(n3)=z(n3+1);end n5=n5+1;endif n1=1&n2=1for n3=1:(n-2) x1(n3)=z(n3+1);end n5=n5+2;endif n 1=0&n 2=1for n3=1:( n-1)x1(n 3)=z( n3);endn5=n5+1;endx=x1;if n 1=0&n 2=0n4=1;enden

24、d第十四周：1.大肠杆菌比生长速率测定。在一定培养条件下，培养大肠杆菌，测得实验数据如下表。求：该条件下，大肠杆菌的最大 比生长速率卩m,半饱和常数Ks，并作模型检验。S (mg/L )卩(h-1 )S (mg/L )u (h-1 )60.061220.60130.121530.66330.241700.69400.312210.70640.432100.731020.53s=6 13 33 40 64 102 122 153 170 221 210;mu=0.06 0.12 0.24 0.31 0.43 0.53 0.60 0.66 0.69 0.70 0.73; spmu=s./mu ;n

25、=len gth(s);a=polyfit(s,spmu,1);mum=1/a(1)ks=a(2)/a(1)lxx=sum(s.A2)-1/ n*(sum(s)f2;lyy=sum(spmu.A2)-1/ n*(sum(spmu)F2;lxy=sum(s.*spmu)-1/ n*sum(s)*sum(spmu); r=lxy/(sqrt(lxx*lyy)R=corrcoef(s,spmu)Qr=lxyA2/lxx;Q=(lxx*lyy-lxyA2)/lxx;F=Qr/(Q/(n-2)2. 多元线性回归Pa=9.0 8.6 8.4 7.5 7.0 6.8 6.5 6.0; Pb=8.3 7.0

26、6.2 4.2 3.9 3.5 2.6 2.2;Pc=2.7 4.4 5.4 8.3 9.1 9.7 10.9 11.8; r=1.97 1.05 0.73 0.25 0.18 0.13 0.07 0.04; k0=ones(8,1);alpha=0.05;r0=log(r);Pa0=log(Pa); Pb0=log(Pb);Pc0=log(Pc); p=k0 Pa0 Pb0 Pc0;b,bint,r,rint,stats=regress(r0,p,alpha)k=exp(b(1) m=r*r p1=Pa0 Pb0 Pc0;stepwise(p1,r0)第十六周1. 对作业( 7 )的两题，分别作非线性回归，并比拟参数值和残差。 function y=ecolinlin(beta0)s=6 13 33 40 64 102 12