南京邮电大学数学实验MATLAB2023

声明

使用R2023b版本完成，参考教材为数学实验（第三版）金正猛 王正新等编 科学出版社。参考了不少前人智慧的结晶，在此表示感谢，不保证全对，仅供大家参考，请勿抄袭，有错漏、疑问之处，烦请通过QQ或者评论区留言交流指正，欢迎大家多交流。联系QQ：3468039783

模块一：基础练习

1.1

clearsyms xm=628;f=(log(1+x-m*x^2)-x)/(1-cos(x));limit(f)f=((sqrt(m*x^2+2)-atan(m*x))/x);limit(f,x,inf)pretty(limit(f,x,inf))

常用函数见书P5页，微积分部分见书P25页。pretty（f）将f显示为数学书写形式，见课本P13页 

1.2

clearsyms xm=628;f=exp(m*x)*sin(x);diff(f,2)subs(f,x,0);diff(f,6)

1.3

clearsyms xm=628;f=(x+sin(x))/(1+cos(x));int(f,x)pretty(int(f,x))f=log(1+m*x)-m*x;int(f,x,0,1)pretty(int(f,x,0,1))

1.4

clearsyms xm=628;f=cos(x*(m/200+sin(x)));taylor(f,x,'Order',5,'ExpansionPoint',0)pretty(taylor(f,x,'Order',5,'ExpansionPoint',0))

在学校的机房的电脑上可能不可以加'Order'，答案可能会变得很长，请自行取舍'Order','ExpansionPoint',0等

1.5

clearx=[];x(1)=rand(1);m=628;for k=2:10    x(k)=sqrt(m/100+x(k-1));endx

rand函数见书P5页

1.6 

clearm=628;A=[4,2,m-2;-3,0,5;1,5,2*m]B=[3,4,0;2,0,-3;-2,1,1]det(A)inv(A)eig(A)[P,D]=eig(A);X=P(1:3,1:1)Y=P(1:3,2:2)Z=P(1:3,3:3)A\B  A/BC=[A B];Q=rref(C)

向量与矩阵见书P11、12页，线性代数见书P27页

1.7

1.7.1

f.m

function y=f(x)if 0<=x&&x<=1/2    y=2*x;elseif 1/2<x&&x<=1    y=2*(1-x);end

g.m

function y=g(x,f)n=length(x);for i=1:n    y(i)=f(x(i));endend

main

fplot(@(x)g(x,@f),[0,1])

建议在自己电脑上做的可以把代码都新建.m文件、规范命名并保存到合适的路径

关于M文件与函数的定义，请翻看书P14页，程序流程控制见书P16页，length函数在书P11页，MATLAB软件中的做题，请从书P30页开始看，注意P33页的常用图形控制命令

保存f.m和g.m两个文件，并运行添加路径，然后在命令行窗口中输入main中的代码，回车执行

1.7.2

f.m

function y=f(x)if -pi<=x&&x<=0    y=x-1;elseif 0<x&&x<=pi    y=x+1;end

g.m

function y=g(x,f)n=length(x);for i=1:n    y(i)=f(x(i));endend

main

fplot(@(x)g(x,@f),[-pi,pi])

保存f.m和g.m两个文件，并运行添加路径，然后在命令行窗口中输入main中的代码，回车执行

1.8

1.8.1

clearm=628;t=-m/100:0.1:m/100;x=(m/20)*cos(t);y=(m/20)*sin(t);z=t;plot3(x,y,z)grid on

注意m值不同时取的t的值也不同

1.8.2

clearm=628;t=-m/100:0.1:m/100;x=cos(t)+t.*sin(t);y=sin(t)-t.*cos(t);z=-t;plot3(x,y,z)grid on

注意m值不同时取的t的值也不同

1.9

clearm=628;x=-30:0.1:30;t=[1000/m,500/m,100/m];for i=1:3    j=t(i);    y=1/(sqrt(2*pi*j))*exp(-(x.^2)/2*(j^2));    plot(x,y)    hold onend

1.10

clearezplot('log(x^2+y*628)-(x^3)*y-sin(x)')

这儿请把m写进式子里

1.11

clearm=628;x=-5:0.1:5;y=-5:0.1:5;[x,y]=meshgrid(x,y);z=m.*x.*exp(-(x.^2+y.^2));mesh(x,y,z)

1.12

clearsyms xm=628;y=x^3+sqrt(m)*x^2+(m/3-3)*x-sqrt(m)*(1-m/27);fplot(x,y,[-sqrt(m)/3-2,-sqrt(m)/3+2])axis([-10.5,-6,-2.5,2.5])grid onfzero('x^3+sqrt(628)*x^2+(628/3-3)*x-sqrt(628)*(1-628/27)',-10)fzero('x^3+sqrt(628)*x^2+(628/3-3)*x-sqrt(628)*(1-628/27)',-8)fzero('x^3+sqrt(628)*x^2+(628/3-3)*x-sqrt(628)*(1-628/27)',-6)diff(y)figurefplot(x,diff(y),[-sqrt(m)/3-2,-sqrt(m)/3+2])axis([-10,-7,-3,1])grid onfzero('3*x^2 + 4*157^(1/2)*x + 619/3',-10)fzero('3*x^2 + 4*157^(1/2)*x + 619/3',-7)

画出图形后看图估计零点的值，然后用fzero求出近似解（fezro里要把m代进去）。再对原函数求导，得出的式子再重复上一步。即可得到单调区间

模块二：函数的迭代

2.1

clearsyms xm=628;y=(2*x+1)/(x-m)-x;solve(y,x)

用以上代码，求解f(x)的两个不动点

编写普适性迭代的MATLAB函数

dd.m

function p=dd(f,x0,n)p=x0;for i=2:n    p(i)=f(p(i-1));endend

main

dd(@(x)(2*x+1)/(x-m),___,20)

在命令行窗口中输入main中的代码，回车执行。在___中填入两个不动点等不同的初值，得出结论，关于超稳定点合吸引点、排斥点，请翻看书P64、65页

2.2

clearm=628;f=@(x)(1-2*abs(x-1/2));t=[1/4,1/10,1/100,1/m];for j=1:4    subplot(2,2,j)    x0=t(j);    for i=1:1:100    plot(i,f(x0),'+');    x0=f(x0);    hold on    endend

可以用F10逐步步进，观察离散点图的绘制过程

2.3

Martin.m

function Martin(a,b,c,N)m=628;f=@(x,y)(y-sign(x)*sqrt(abs(b*x-c)));g=@(x)(a-x);m=[0;0];for n=1:N    m(:,n+1)=[f(m(1,n),m(2,n)),g(m(1,n))];endplot(m(1,:),m(2,:),'kx')axis equal

main

m=628subplot(2,2,1)Martin(m,m,m,5000)subplot(2,2,2)Martin(m/1000,m/1000,0.5,10000)subplot(2,2,3)Martin(m/1000,m,-m,15000)subplot(2,2,4)Martin(-m/10,17,4,20000)

在命令行窗口中输入main中的代码，回车执行

2.4

clearm=628;syms x a b c d e;diff((a*x+m*c)/(c*x^2+a))pretty(diff((a*x+m*c)/(c*x^2+a)))x=m^(1/3)vpa(diff(subs((100*x+m)/(x^2+100),x,m^(1/3))),8) f=@(x)(100*x+628)/(x^2+100);for j=1:5    x0=(100*rand-50)    for i=1:20        x0=f(x0);        fprintf('%g,%g\n',i,x0)    end end

这儿abcde的值需要先将f(x)=x求解，得出a=e,d=0,c*x^3=b,又x=m^(1/3),故b=cm,原式化简为(a*x+m*c)/(c*x^2+a)。此处，我假设a=100，c=1。x=m^(1/3)约等于8，用diff求解导数，在x=m^(1/3)时导数小于1，满足条件。(-50,50)范围内取5个随机数，迭代20次，观察x0的变化情况。可以更改x0=(100*rand-50)尝试在不同数量级下迭代的情况

2.5

clearsyms xy=sin(x);y1=taylor(y,x,'Order',2,'ExpansionPoint',0);y2=taylor(y,x,'Order',4,'ExpansionPoint',0);y3=taylor(y,x,'Order',6,'ExpansionPoint',0);y4=taylor(y,x,'Order',8,'ExpansionPoint',0);y5=taylor(y,x,'Order',10,'ExpansionPoint',0);y6=taylor(y,x,'Order',12,'ExpansionPoint',0);subplot(1,2,1)fplot([y,y1,y2,y3],[-3*pi/2,3*pi/2])grid onsubplot(1,2,2)fplot([y,y4,y5,y6],[-3*pi/2,3*pi/2])grid on

结论请参见书P99页

模块三:线性映射的迭代

3.1、3.2

clearA=str2sym('[628,628-4;6-628,10-628]');[P,D]=eig(A);Q=inv(P);syms n;x=[1;2];xn=P*(D.^n)*Q*xB=1/10.*A[P,D]=eig(B);Q=inv(P);xn=P*(D.^n)*Q*x

得把m的值代进去，不然结果中还会有m

3.3

clearA=[9,5;2,6]for k=1:20    x=[2*rand-1;2*rand-1]    t=[];    for i=1:40        x=A*x;        t(i,1:2)=x;    end    for j=1:40        if t(j,1)==0        else a=t(j,2)/t(j,1);            fprintf('%g,%g\n',j,a);        end    end    plot(t(1:40,1),t(1:40,2),'+')    grid on    hold onend

结论请参见书P136页

3.4

clearm=628;A=[3/4,7/18;1/4,11/18];A(:,:,2)=[m,6-m;m-2,8-m];A(:,:,3)=[m,1/4-m;m-3/4,1-m];A(:,:,4)=[m-1,m;1-m,-m];for i=1:4eig(A(:,:,i))[P,D]=eig(A(:,:,i))X=P(1:2,1:1)Y=P(1:2,2:2)    for j=1:20        x=[2*rand-1;2*rand-1];        t=[];        for k=1:40            x=A(:,:,i)*x;            t(k,1:2)=x;        end        t        subplot(2,2,i)        plot(t(1:40,1),t(1:40,2),'+')        grid on        hold on    endend

这儿我参照了3.3，对每个迭代矩阵都取了20组随机数作为初始向量，实际上只需要取x=[0.5,0.5]即可，请自行修改代码，对于迭代后趋向于无穷大(零)的，可能需要归一化（见书P138页）

模块四：种群增长模型与综合实验

4.1、4.2

clearf=[1,2,4,5]for i=1:4    for b=1:1000        a=sqrt((b+f(i))^2-b^2);        if (a==floor(a))            fprintf(['a=%4i,b=%4i,c=%4i|'],a,b,b+f(i))        end    endfprintf('\n')end

详见书P159-162页

4.3

clearfor k=1:200 for b=1:1000a=sqrt((b+k)^2-b^2);if((a==floor(a))&&gcd(gcd(a,b),(b+k))==1)            fprintf('%i,',k);break;endendend

4.4

clearsyms a0 a1m=9;x=[1.5,1.8,2.4,2.8,3.4,3.7,4.2,4.7,5.3];y=[8.9,10.1,12.4,14.3,16.2,17.8,19.6,22.0,24.1];polyfit(x,y,1)x1=sum(x);x2=sum(x.^2);y1=sum(y);xy=sum(y.*x);eql1=m*a0+x1*a1==y1eql2=x1*a0+x2*a1==xy[a0,a1]=solve([eql1,eql2],[a0,a1])error=sum((y-(a0+a1.*x)).^2)f=a0+a1*x;plot(x,y,'+',x,f);

关于最小二乘法和法方程组，请参见书P192、193页

4.5

clearsyms k x0X=1790:10:1980;x=1790:10:1930;Y=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204.0,226.5];y=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,106.5,123.2];Z=log(Y);z=log(y);A=[1,x(1);1,x(15)];B=[z(1);z(15)];u=A\B;x0=u(1)k=u(2)error=sum((x0*exp(k*x)-y).^2)f=exp(x0+k*X);plot(X,f,X,Y,'+');hold ont=polyfit(X,Z,1);X0=t(2)K=t(1)f=exp(X0+K*X);error=sum((X0*exp(K*X)-Y).^2)plot(X,f);

先使用了两个点求解，再使用所有的数据用最小二乘的方法求解，error这一式子不一定正确

4.6

clearx=1:26;y=[1807,2001,2158,2305,2422,2601,2753,2914,3106,3303,3460,3638,3799,3971,4125,4280,4409,4560,4698,4805,4884,4948,5013,5086,5124,5163];a=[6000,2,0.1];f=@(a,x)a(1)./(1+a(2)*exp(-a(3)*x));[A,resnorm]=lsqcurvefit(f,a,x,y);t=26;while abs(f(A,t)-f(A,t+1))>1    t=t+1;endf(A,t)tt=1:t;plot(x,y,'+',t,f(A,t))

4.7

clearx=1:26;y=[1807,2001,2158,2305,2422,2601,2753,2914,3106,3303,3460,3638,3799,3971,4125,4280,4409,4560,4698,4805,4884,4948,5013,5086,5124,5163];a=[6000,2,0.1,0.1];f=@(a,x)a(1)./(1+a(2)*exp(-a(3)*x-a(4)*x.^2));[A,resnorm]=lsqcurvefit(f,a,x,y);t=26;while abs(f(A,t)-f(A,t+1))>1    t=t+1;endf(A,t)tt=1:t;plot(x,y,'+',t,f(A,t));

4.8

clearx=1:27;y=[21,65,127,281,558,923,1321,1801,2420,3125,3886,4638,5306,6028,6916,7646,8353,9049,9503,10098,10540,10910,11287,11598,11865,12086,12251];a=[13000,2,0.1];f=@(a,x)a(1)./(1+a(2)*exp(-a(3)*x));[A,resnorm]=lsqcurvefit(f,a,x,y);t=27;while abs(f(A,t)-f(A,t+1))>1    t=t+1;endf(A,t)tt=1:t;plot(x,y,'+',t,f(A,t));hold ona=[13000,2,0.1,0.1];f=@(a,x)a(1)./(1+a(2)*exp(-a(3)*x-a(4)*x.^2));[A,resnorm]=lsqcurvefit(f,a,x,y);t=27;while abs(f(A,t)-f(A,t+1))>10    t=t+1;endf(A,t)tt=1:t;plot(t,f(A,t));

使用4.7中的模型时，while语句中的判断范围不能是1，否则模型有误，可以更改为10

4.9

clearx=1:27;y=[21,65,127,281,558,923,1321,1801,2420,3125,3886,4638,5306,6028,6916,7646,8353,9049,9503,10098,10540,10910,11287,11598,11865,12086,12251];a=[2,0.1,0.1];f=@(a,x)a(1)*exp(exp(a(2)*x+a(3)));[A,resnorm]=lsqcurvefit(f,a,x,y);t=27;while abs(f(A,t)-f(A,t+1))<1    t=t+1;endf(A,t)tt=1:t;plot(x,y,'+',t,f(A,t));

综合实验

clearfor j=1:25    p=rand;    AA=p*p;    Aa=2*p*(1-p);    aa=(1-p)*(1-p);    K=AA+Aa+aa;    for i=1:2        x(i)=AA(i)*AA(i)+AA(i)*Aa(i)*1/2+Aa(i)*AA(i)*1/2+Aa(i)*Aa(i)*1/4;        y(i)=AA(i)*Aa(i)*1/2+AA(i)*aa(i)+Aa(i)*AA(i)*1/2+Aa(i)*Aa(i)*1/2+Aa(i)*aa(i)*1/2+aa(i)*AA(i)+aa(i)*Aa(i)*1/2;        z(i)=aa(i)*aa(i)+Aa(i)*Aa(i)*1/4+Aa(i)*aa(i)*1/2+aa(i)*Aa(i)*1/2;        k=x+y+z;        AA(i+1)=x(i);        Aa(i+1)=y(i);        aa(i+1)=z(i);    end    tu=[AA,;Aa,;aa,];    subplot(5,5,j);    plot(1:3,tu');    grid on    title(j);    legend('AA','Aa','aa');end

我选择的是第二个实验，这是我最终版的代码，前两个版本的代码在压缩包里。