matlab數(shù)學(xué)建模實(shí)例

1、第四周3.function y=mj()for x0=0:0.01:8 x1=x03-11.1*x02+38.79*x0-41.769; if (abs(x1)<1.0e-8) x0 endend4. 分別用簡單迭代法、埃特金法、牛頓法求解方程，并比較收斂性與收斂速度（e分別取10-3、10-5、10-8）。簡單迭代法：function y=jddd(x0)x1=(20+10*x0-2*x02-x03)/20; k=1;while (abs(x1-x0)>=1.0e-3) x0=x1; x1=(20+10*x0-2*x02-x03)/20;k=k+1;endx1k埃特金法：func

2、tion y=etj(x0)x1=(20-2*x02-x03)/10;x2=(20-2*x12-x13)/10;x3=x2-(x2-x1)2/(x2-2*x1+x0);k=1;while (abs(x3-x0)>=1.0e-3) x0=x3; x1=(20-2*x02-x03)/10; x2=(20-2*x12-x13)/10; x3=x2-(x2-x1)2/(x2-2*x1+x0);k=k+1;endx3k牛頓法：function y=newton(x0)x1=x0-fc(x0)/df(x0);k=1;while (abs(x1-x0)>=1.0e-3) x0=x1; x1=x0

3、-fc(x0)/df(x0);k=k+1;endx1kfunction y=fc(x)y=x3+2*x2+10*x-20;function y=df(x)y=3*x2+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*x

4、0+G;n=1;while norm(s-x0)>=1.0e-8 x0=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*U;G=C*d;s=B*x0+G;n=1;while norm(s-x0)>=1.0e-5 x0=s; s=B*x0+G; n=n+1;endn松弛法：function s=loose(a,d,x0,w)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1)

5、;C=inv(D-w*L);B=C*(1-w)*D+w*U);G=w*C*d;s=B*x0+G;n=1;while norm(s-x0)>=1.0e-8 x0=s; s=B*x0+G; n=n+1;endn2.練習(xí)MATLAB的常用矩陣語句，就龍格現(xiàn)象函數(shù)（p88）練習(xí)插值語句interp, spline，并比較。 3.測得血液中某藥物濃度隨時(shí)間的變化值為: 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 時(shí)的濃

6、度C. 分別用n=4,5,9的拉格朗日插值計(jì)算；并用樣條函數(shù)插值計(jì)算,并比較結(jié)果。拉格朗日插值:function s=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.88;x0=0.45 1.75 5.0 6.0;m=length(x0);for i=1:m D=abs(x-x0(i); I=1; while I<=n+1 for a=1:length(x) if D(a)=min(D) c(I)=a; D(a)=max(D)+

7、1; break end end I=I+1; end b=sort(c); z=x0(i); t=0.0; for k=1:length(b) u=1.0; for j=1:length(b) if j=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樣條函數(shù)差值：Interp1(x,y,x0,spline)Spline(x,y,x0)第八周1.給定某藥物濃度隨時(shí)間的變化值(作業(yè)3)，1）分別采用樣條函數(shù)和三點(diǎn)公式（設(shè)h=0.1)求結(jié)點(diǎn)處的導(dǎo)數(shù)值，并比較結(jié)果。2）求該時(shí)間段的平均濃度（定步長S法）樣條函數(shù)

8、：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)三點(diǎn)公式：function df=sandian()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;h=0.1;n=

9、length(t);for i=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); switch i case 1 df(i)=(-3*y0+4*y1-y2)/(2*h); case n df(i)=(y4-4*y3+3*y0)/(2*h); otherwise df(i)=(y1-y3)/(2*h); endendend平均濃度：function averagec=simpson()t=0.25 0.5 1.0 1.5 2.0

10、 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.練習(xí)MATLAB常用的trapz, quad, quadl等語句。計(jì)算：x=0:8;y=1./(sqrt(2.*pi).*exp(-(x-4).2./2);z=trapz(x,y)function y=jifen(x)y=1./(sqrt(2.*pi).*exp(-(x-4).2./2);q1=quad('

11、;jifen',0,8,1.0e-8)q2=quadl('jifen',0,8,1.0e-8)3.采用變步長經(jīng)典R-K法, ode23, ode45計(jì)算例9-5，并作比較。變步長經(jīng)典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=

12、1:m1 k1=prey(t,x); for i=1:n x2(i)=x(i)+1/2*h*k1(i); end k2=prey(t+h/2,x2); for i=1:n x3(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); end k4=prey(t+h,x4); for i=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;endh1=length(z);t2=a:(b-a)/(h1-1):b;plo

13、t(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=prey(t0+h/4,x0+h/4*k12);k32=prey(t0+h/4,x0+h/4*k22);k

14、42=prey(t0+h/2,x0+h/2*k32);x2=x0+h/12*(k12+2*k22+2*k32+k42);if abs(x2-x1)<1.0e-5 while abs(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

15、=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 while abs(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

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

17、)t,x=ode45('prey',0:0.1:15,25 2);plot(t,x)t,xgtext('x1(t)')gtext('x2(t)')第十周1熟悉常用的概率分布、概率密度函數(shù)圖、分位點(diǎn)。（統(tǒng)計(jì)工具箱）2對例10-1作統(tǒng)計(jì)分組(每組間隔分別為3cm、5cm)，并作直方圖，計(jì)算特征值與置信區(qū)間；如假設(shè) 0=175作檢驗(yàn)（=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 182 160 159

18、172 178 166 159 173 161 150 164 175 173 163 165 146 163 162 158 164 169 170 164 179 169 178 170 155 169 160 174 159 168 151 176 164 161 163 172 167 154 164 153 165 161 168 166 166 148 161 163 177 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

19、179 165 184 157;m=ceil(max(data)-min(data)/n);hist(data,m)data=162 166 171 167 157 168 164 178 170 152 158 153 160 174 159 167 171 168 182 160 159 172 178 166 159 173 161 150 164 175 173 163 165 146 163 162 158 164 169 170 164 179 169 178 170 155 169 160 174 159 168 151 176 164 161 163 172 167 154 1

20、64 153 165 161 168 166 166 148 161 163 177 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自行尋找生物學(xué)數(shù)據(jù)，進(jìn)行分析，試作曲線圖、條形圖、餅圖。（包括圖示） 第十二周1、作圖練

21、習(xí)不同形式誤差的疊加，隨機(jī)誤差+周期性誤差；隨機(jī)誤差+線性誤差；隨機(jī)誤差+恒定誤差。(自行設(shè)計(jì)數(shù)據(jù)，注意誤差數(shù)量級的選?。?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('

22、;S/(g.kg-1 of water)')title('Solubility of ¦Á-Form Glycine in Water') 3、異常數(shù)據(jù)剔除 拉依特準(zhǔn)則：function y=lyt()x=25.307 25.112 25.324 25.300 25.295 25.293 25.294 25.314 25.341 25.315 25.314 25.299 25.303 25.313 25.311 25.590 25.309 25.316 25.310 25.317 25.306 25.291 25.325 25.010 25.315

23、25.438;mu=mean(x);sigma=std(x);n=length(x);if n<10 m=2;else m=3;for i=1:n if abs(x(i)-mu)>m*sigma i x(i) endendend格魯布斯準(zhǔn)則：function y=grubbs()x=25.307 25.112 25.324 25.300 25.295 25.293 25.294 25.314 25.341 25.315 25.314 25.299 25.303 25.313 25.311 25.590 25.309 25.316 25.310 25.317 25.306 25.29

24、1 25.325 25.010 25.315 25.438;mu=mean(x);sigma=std(x);n=length(x);for i=1:n if abs(x(i)-mu)/sigma)>=2.681 %格魯布斯極限值（n=26）：0.005-3.157 0.01-3.029 0.025-2.841 0.05-2.681 i x(i) endendend狄克遜準(zhǔn)則：x=25.307 25.112 25.324 25.300 25.295 25.293 25.294 25.314 25.341 25.315 25.314 25.299 25.303 25.313 25.311 2

25、5.590 25.309 25.316 25.310 25.317 25.306 25.291 25.325 25.010 25.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=0 z=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) end n2=0; b=(z(n)-z(n-2)/(z(n)-z(3); if

26、b>f(n5) n2=1; z(n) end x1=0 0; if n1=1&&n2=0 for n3=1:(n-1) x1(n3)=z(n3+1); end n5=n5+1; end if n1=1&&n2=1 for n3=1:(n-2) x1(n3)=z(n3+1); end n5=n5+2; end if n1=0&&n2=1 for n3=1:(n-1) x1(n3)=z(n3); end n5=n5+1; end x=x1; if n1=0&&n2=0 n4=1; endend第十四周：1.大腸桿菌比生長速率測定

27、。在一定培養(yǎng)條件下，培養(yǎng)大腸桿菌，測得實(shí)驗(yàn)數(shù)據(jù)如下表。求：該條件下，大腸桿菌的最大比生長速率m，半飽和常數(shù)Ks，并作模型檢驗(yàn)。 S（mg/L） （h-1） S（mg/L） （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=length(s);a=polyfit(s,spmu

28、,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=lxy2/lxx;Q=(lxx*lyy-lxy2)/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 6.2 4.2 3.9 3.5 2.6 2.2'Pc

29、=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'*rp1=Pa0 Pb0 Pc0;stepwise(p1,r0)第十六周1. 對作業(yè)（7）的兩題，分別作非線性回歸，并比較參數(shù)值和殘差。 function y=ecolinlin(beta0)s=6 13 33