matlab数学建模实例

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-3) x0=x1; x1=(20+10*x0-2*x02-x03)/20;k=k+1;endx1k埃特金法：function 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-

2、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-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

3、=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-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

4、+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);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.练习MATLAB的常用矩阵语句，就龙格现象函数（p88）练习插值语句interp, spline，并比较。 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的拉格朗日插值计算；并用样条函数插值计算,并比较结果。拉格朗日插值: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

6、.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)+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样条函数差值：In

7、terp1(x,y,x0,spline)Spline(x,y,x0)第八周1.给定某药物浓度随时间的变化值(作业3)，1）分别采用样条函数和三点公式（设h=0.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(

8、)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=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*

9、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 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./

10、(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(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=ze

11、ros(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: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)

12、+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;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,x

13、0);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);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+k4

14、1); 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); 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

15、)t,x=ode45(prey,0:0.1:15,25 2);plot(t,x)t,xgtext(x1(t)gtext(x2(t)第十周1熟悉常用的概率分布、概率密度函数图、分位点。（统计工具箱）2对例10-1作统计分组(每组间隔分别为3cm、5cm)，并作直方图，计算特征值与置信区间；如假设 0=175作检验（=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 172 178 166 159 173 161 150 164

16、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 179 165 184 157;m=ceil(max(data)

17、-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 164 153 165 161 168 166 166 148 1

18、61 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自行寻找生物学数据，进行分析，试作曲线图、条形图、饼图。（包括图示） 第十二周1、作图练习不同形式误差的叠加，随机误差+周期性误差；随机误差+线性误差；

19、随机误差+恒定误差。(自行设计数据，注意误差数量级的选取） 2、作errorbar图（本课件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 -Form Glycine in Water) 3、

20、异常数据剔除 拉依特准则：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 25.438;mu=mean(x);sigma=std(x);n=length(x);if nm*sigma i x(i) endendend格鲁布斯准则：function y=grubb

21、s()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 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

22、0.025-2.841 0.05-2.681 i x(i) endendend狄克逊准则：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 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.50

23、7 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 af(n5) n1=1; z(1) end n2=0; b=(z(n)-z(n-2)/(z(n)-z(3); if bf(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

24、 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.大肠杆菌比生长速率测定。在一定培养条件下，培养大肠杆菌，测得实验数据如下表。求：该条件下，大肠杆菌的最大比生长速率m，半饱和常数Ks，并作模型检验。 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

25、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,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-lxy

26、2)/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=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,al

27、pha)k=exp(b(1)m=r*rp1=Pa0 Pb0 Pc0;stepwise(p1,r0)第十六周1. 对作业（7）的两题，分别作非线性回归，并比较参数值和残差。 function y=ecolinlin(beta0)s=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;beta,r,J=nlinfit(s,mu,szsl,beta0);mum=beta(1)ks=beta(2)rJfunction mu=szsl(beta,s)mu=beta(

28、1).*s./(beta(2)+s);endfunction y=reacnlin(beta0)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=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;p=Pa Pb Pc;beta,r,J=nlinfit(p,R,fydl,beta0);k=beta(4)n1=beta(1)n2=beta(2)n3=beta(3)rJfunction r=fydl(be

29、ta,p)r=beta(4).*p(:,1).beta(1).*p(:,2).beta(2).*p(:,3).beta(3);end2. 作二次正交回归。数据如右，比较不同模型计算结果。x1=1 1 1 1 -1 -1 -1 -1 -1.68 1.68 0 0 0 0 0 0 0 0 0 0;x2=1 1 -1 -1 1 1 -1 -1 0 0 -1.68 1.68 0 0 0 0 0 0 0 0;x3=1 -1 1 -1 1 -1 1 -1 0 0 0 0 -1.68 1.68 0 0 0 0 0 0;y=730.2 780.5 266.7 224.5 783.1 837.5 622.6 5

30、38.3 536.2 221.2 214.2 926.2 702.4 680.1 868.5 788.3 856.5 853.4 772.6 848.4;x=x1 x2 x3;alpha=0.05;rstool(x,y,linear,alpha)附录资料：MATLAB的30个方法1 内部常数pi 圆周率 exp(1)自然对数的底数ei 或j 虚数单位Inf或 inf 无穷大 2 数学运算符a+b 加法a-b减法a*b矩阵乘法a.*b数组乘法a/b矩阵右除ab矩阵左除a./b数组右除a.b数组左除ab 矩阵乘方a.b数组乘方-a负号 共轭转置.一般转置3 关系运算符=等于大于=大于或等于=不等于

31、4 常用内部数学函数 指数函数exp(x)以e为底数对数函数log(x)自然对数，即以e为底数的对数log10(x)常用对数，即以10为底数的对数log2(x)以2为底数的x的对数开方函数sqrt(x)表示x的算术平方根绝对值函数abs(x)表示实数的绝对值以及复数的模三角函数（自变量的单位为弧度）sin(x)正弦函数cos(x)余弦函数tan(x)正切函数cot(x)余切函数sec(x)正割函数csc(x)余割函数反三角函数 asin(x)反正弦函数acos(x)反余弦函数atan(x)反正切函数acot(x)反余切函数asec(x)反正割函数acsc(x)反余割函数双曲函数 sinh(x)

32、双曲正弦函数cosh(x)双曲余弦函数tanh(x)双曲正切函数coth(x)双曲余切函数sech(x)双曲正割函数csch(x)双曲余割函数反双曲函数 asinh(x)反双曲正弦函数acosh(x)反双曲余弦函数atanh(x)反双曲正切函数acoth(x)反双曲余切函数asech(x)反双曲正割函数acsch(x)反双曲余割函数求角度函数atan2(y,x)以坐标原点为顶点，x轴正半轴为始边，从原点到点（x，y）的射线为终边的角，其单位为弧度，范围为（ ， 数论函数gcd(a,b)两个整数的最大公约数lcm(a,b)两个整数的最小公倍数排列组合函数factorial(n)阶乘函数，表示n的

33、阶乘 复数函数 real(z)实部函数imag(z)虚部函数abs(z)求复数z的模angle(z)求复数z的辐角，其范围是（ ， conj(z)求复数z的共轭复数求整函数与截尾函数ceil(x)表示大于或等于实数x的最小整数floor(x)表示小于或等于实数x的最大整数round(x)最接近x的整数最大、最小函数max(a，b，c，)求最大数min(a，b，c，)求最小数符号函数 sign(x)5 自定义函数-调用时：“返回值列=M文件名（参数列）”function 返回变量=函数名（输入变量） 注释说明语句段（此部分可有可无）函数体语句 6进行函数的复合运算compose(f,g) 返回值

34、为f(g(y)compose(f,g,z) 返回值为f(g(z)compose(f,g,x,.z) 返回值为f(g(z)compose(f,g,x,y,z) 返回值为f(g(z)7 因式分解syms 表达式中包含的变量 factor(表达式) 8 代数式展开syms 表达式中包含的变量 expand(表达式)9 合并同类项syms 表达式中包含的变量 collect(表达式，指定的变量)10 进行数学式化简syms 表达式中包含的变量 simplify(表达式)11 进行变量替换syms 表达式和代换式中包含的所有变量 subs(表达式，要替换的变量或式子，代换式)12 进行数学式的转换调用M

35、aple中数学式的转换命令，调用格式如下：maple(Maple的数学式转换命令) 即：maple(convert(表达式，form)将表达式转换成form的表示方式 maple(convert(表达式，form, x) 指定变量为x，将依赖于变量x的函数转换成form的表示方式（此指令仅对form为exp与sincos的转换式有用） 13 解方程solve(方程，变元) 注：方程的等号用普通的等号： = 14 解不等式调用maple中解不等式的命令即可，调用形式如下： maple(maple中解不等式的命令)具体说，包括以下五种：maple( solve（不等式）) maple( solve

36、（不等式，变元） ) maple( solve（不等式，变元） ) maple( solve（不等式，变元） ) maple( solve（不等式，变元） )15 解不等式组调用maple中解不等式组的命令即可，调用形式如下： maple(maple中解不等式组的命令) 即：maple( solve（不等式组，变元组） )16 画图方法：先产生横坐标的取值和相应的纵坐标的取值，然后执行命令： plot(x,y) 方法2：fplot(f(x)，xmin,xmax) fplot(f(x)，xmin,xmax,ymin,ymax) 方法3：ezplot(f(x) ezplot(f(x) ,xmin,

37、xmax) ezplot(f(x) ,xmin,xmax,ymin,ymax) 17 求极限（1）极限：syms x limit(f(x), x, a) （2）单侧极限：左极限：syms x limit(f(x), x, a,left)右极限：syms x limit(f(x), x, a,right) 18 求导数diff(f(x) diff(f(x),x) 或者：syms x diff(f(x) syms x diff(f(x), x) 19 求高阶导数 diff(f(x)，n) diff(f(x),x，n)或者：syms x diff(f(x)，n)syms x diff(f(x), x，n) 20 在MATLAB中没