maab数学建模实例

1、3 .function y=mj()for x0=0:8x1=x0A*x0A2+*;if (abs(x1)<x0endend4 .分别用简单迭代法、埃特金法、牛顿法求解方程，并比较收敛性与收敛速度(分别取 10-3、10-5、10-8)。简单迭代法：function y=jddd(x0)x1=(20+10*x0-2*x0A2-x0A3)/20;k=1;while (abs(x1-x0)>=x0=x1;x1=(20+10*x0-2*x0A2-x0A3)/20;k=k+1;endx1k埃特金法：function y=etj(x0)x1=(20-2*x0A2-x0A3)/10;x2=(2

2、0-2*x1A2-x1A3)/10;x3=x2-(x2-x1)A2/(x2-2*x1+x0);k=1;while (abs(x3-x0)>=x0=x3;x1=(20-2*x0A2-x0A3)/10;x2=(20-2*x1A2-x1A3)/10;x3=x2-(x2-x1)A2/(x2-2*x1+x0);k=k+1;endx3k牛顿法:function y=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 kfunction y=fc(x) y=xA

3、3+2*xA2+10*x-20; function y=df(x) y=3*xA2+4*x+10; 第六周1 .解例6-4 (p77)的方程组，分别采用消去法(矩阵分解)、Jacobi迭代法、Seidel迭代法、松弛法求解，并比较收敛速度。消去法： 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)>= x0=s

4、;s=B*x0+G;n=n+1; end nSeidel迭代法：function s=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;while norm(s-x0)>=x0=s;s=B*x0+G;n=n+1;end n松弛法：function s=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

5、;n=1;while norm(s-x0)>=x0=s;s=B*x0+G;n=n+1;end n2 .练习MATLA的常用 矩阵语句，就龙格现象函数(p88)练习插值语句interp, spline , 并比较。分别用n=4,5,9的拉格朗日插值计算；并用样条函数插值计算，并比较结果 拉格朗日插值:functions=lagr(n)x=；y=；x0=； m=length(x0);for i=1:mD=abs(x-x0(i)；I=1;while I<=n+1for a=1:length(x)if D(a)=min(D) c(I)=a； D(a)=max(D)+1； breakend

6、endI=I+1； endb=sort(c)；z=x0(i)；t=；for k=1:length(b)u=；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,x0, ' spline ')Spline(x,y,x0)第八周卜=求1.给定某药物浓度随时间的变化值(作业3) , 1)分别采用样条函数和三点公式 结点处的导数值，并比较结果。2)求该时间段的平均浓度(定步长St)样条函数：x=;y=;pp=csape(x,y,

7、 'not-a-knot' );df=fnder(pp);df1=ppval(df,x)三点公式：function df=sandian()t=；c=；h=;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); otherwisedf(i)

8、=(y1-y3)/(2*h); endend end 平均浓度： function averagec=simpson() t=;c=;m=(t(1)+t(10)/2;y=spline(t,c,m);averagec=(c(1)+4*y+c(10)/6; end2 .练习MATLA雷用的trapz, quad, quadl 等语句。计算：x=0:8;y=1./(sqrt(2.*pi).*exp(-(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&

9、#39; ,0,8, q2=quadl( 'jifen' ,0,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);k4=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,

10、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);endk3=prey(t+h/2,x3);for i=1:nx4(i)=x(i)+h*k3(i);endk4=prey(t+h,x4);for i=1:nx(i)=x(i)+h/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i);z(i,k)=x(i);endt=t+h;endh1=length(z);t2=a:(b-a)/(h1-1):b;plot(t2,z)gtext( 'x1(t)')gt

11、ext( '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=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=

12、x0+h/12*(k12+2*k22+2*k32+k42);if abs(x2-x1)<while abs(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)

13、;x2=x0+h/12*(k12+2*k22+2*k32+k42);endh=h/2;elsewhile abs(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

14、);x2=x0+h/12*(k12+2*k22+2*k32+k42);endendfunction xdot=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,25 2);plot(t,x)t,xgtext( 'x1(t)')gtext( 'x2(t)')t,x=ode45( 'prey' ,0:15,25 2);plot(t,x)t,xgtext( 'x1(t)')gtext( 

15、9;x2(t)')第十周1 .熟悉常用的概率分布、概率密度函数图、分位点。(统计工具箱)2 .对例10-1作统计分组(每组间隔分别为3cm 5cm),并作直方图，计算特征值与置信区 问；如假设(10=175(乍检验(a =)function y=zf(n)data=162 166 171 167 157 168 164 178 170 152 158 153 160 174 159 167 171 168 182 160 159 172178 166 159 173 161150164175173163 165 146 163 162 158164169170164179169 178

16、 170 155169 160 174 159 168151176164161163 172 167 154 164 153165161168166166148 161 163 177178 171 162 156 165176170156172163 165 149 176 170 182159164179162151170 160 165 167155 168 179 165 184157;m=ceil(max(data)-min(data)/n);hist(data,m)data=162 166 171 167 157 168 164 178 170 152 158 153 160 17

17、4 159 167 171 168 182 160 159 172178 166 159 173 161150164175173163 165 146 163 162 158164169170164179169 178 170 155169 160 174 159 168151176164161163 172 167 154 164 153165161168166166148 161 163 177178 171 162 156 165176170156172163 165 149 176 170 182159164179162151170 160 165 167155 168 179 165

18、 184157;E=mean(data)D=var(data)mu sigma muci sigmaci=normfit(data,h,p,ci=ttest(data,175,0)3 .自行寻找生物学数据，进行分析，试作曲线图、条形图、饼图。(包括图示)第十二周1、作图练习不同形式误差的叠加，随机误差+周期性误差；随机误差+线性误差；随机误差十叵定误差。(自行设计数据，注意误差数量级的选取)2、作 errorbar 图(本课件 Page 3-A)T=;S=;E=;errorbar(T,S,E)xlabel( 'T/')ylabel( 'S/ of water)'

19、 )title( 'Solubility of |& -Form Glycine in Water' )3、异常数据剔除拉依特准则：function y=lyt()x=;mu=mean(x);sigma=std(x);n=length(x);if n<10m=2;else m=3;for i=1:nif abs(x(i)-mu)>m*sigmaix(i)end endend格鲁布斯准则：function y=grubbs()x=；mu=mean(x);sigma=std(x);n=length(x);n=26):for i=1:nif abs(x(i)-mu

20、)/sigma)>=端鲁布斯极限值(ix(i)endendend狄克逊准则：x=;n4=0;f=;while n4=0z=sort(x);n=length(x);n5=1;a=(z(3)-z(1)/(z(n-2)-z(1);n1=0;if a>f(n5)n1=1;z(1)endn2=0;b=(z(n)-z(n-2)/(z(n)-z(3);if b>f(n5)n2=1;z(n)endx1=0 0;if n1=1&&n2=0for n3=1:(n-1)x1(n3)=z(n3+1);endn5=n5+1;endif n1=1&&n2=1for n3=

21、1:(n-2)x1(n3)=z(n3+1);endn5=n5+2;endif n1=0&&n2=1for n3=1:(n-1)x1(n3)=z(n3);endn5=n5+1;endx=x1;if n1=0&&n2=0n4=1;endend第十四周：1 .大肠杆菌比生长速率测定。在一定培养条件下，培养大肠杆菌，测得实验数据如下表。求：该条件下，大肠杆菌的 最大比生长速率以m,半饱和常数Ks,并作模型检验。S (mg/L)/ -1、i (h )S (mg/L)/ -1、i (h )612213153331704022164210102s=6 13 33 40 64

22、102 122 153 170 221 210;mu= ;spmu=s./mu;n=length(s);a=polyfit(s,spmu,1);mum=1/a(1)ks=a(2)/a(1)lxx=sum(s.A2)-1/n*(sum(s)A2;lyy=sum(spmu.A2)-1/n*(sum(spmu)A2;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=,;Pb=,;Pc=,;r=1;k0=ones(8,1);alpha=;r0=log(r);Pa0=log(Pa);Pb0=log(Pb);Pc0=log(Pc);p=k0 Pa0 Pb0 Pc0;b,bint,r,rint,stats=regress(r0