声明

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

模块一：基础练习

1.1

<code>clear

syms x

m=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

clear

syms x

m=628;

f=exp(m*x)*sin(x);

diff(f,2)

subs(f,x,0);

diff(f,6)

1.3

clear

syms x

m=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

clear

syms x

m=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

clear

x=[];

x(1)=rand(1);

m=628;

for k=2:10

x(k)=sqrt(m/100+x(k-1));

end

x

rand函数见书P5页

1.6 

clear

m=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/B

C=[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));

end

end

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));

end

end

main

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

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

1.8

1.8.1

clear

m=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

clear

m=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

clear

m=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 on

end

1.10

clear

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

这儿请把m写进式子里

1.11

clear

m=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

clear

syms x

m=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 on

fzero('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)

figure

fplot(x,diff(y),[-sqrt(m)/3-2,-sqrt(m)/3+2])

axis([-10,-7,-3,1])

grid on

fzero('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

clear

syms x

m=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));

end

end

main

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

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

2.2

clear

m=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

end

end

可以用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))];

end

plot(m(1,:),m(2,:),'kx')

axis equal

main

m=628

subplot(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

clear

m=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

clear

syms x

y=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 on

subplot(1,2,2)

fplot([y,y4,y5,y6],[-3*pi/2,3*pi/2])

grid on

结论请参见书P99页

模块三:线性映射的迭代

3.1、3.2

clear

A=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*x

B=1/10.*A

[P,D]=eig(B);

Q=inv(P);

xn=P*(D.^n)*Q*x

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

3.3

clear

A=[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 on

end

结论请参见书P136页

3.4

clear

m=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:4

eig(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

end

end

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

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

4.1、4.2

clear

f=[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

end

fprintf('\n')

end

详见书P159-162页

4.3

clear

for k=1:200

for b=1:1000

a=sqrt((b+k)^2-b^2);

if((a==floor(a))&&gcd(gcd(a,b),(b+k))==1)

fprintf('%i,',k);

break;

end

end

end

4.4

clear

syms a0 a1

m=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==y1

eql2=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

clear

syms k x0

X=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 on

t=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

clear

x=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;

end

f(A,t)

t

t=1:t;

plot(x,y,'+',t,f(A,t))

4.7

clear

x=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;

end

f(A,t)

t

t=1:t;

plot(x,y,'+',t,f(A,t));

4.8

clear

x=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;

end

f(A,t)

t

t=1:t;

plot(x,y,'+',t,f(A,t));

hold on

a=[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;

end

f(A,t)

t

t=1:t;

plot(t,f(A,t));

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

4.9

clear

x=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;

end

f(A,t)

t

t=1:t;

plot(x,y,'+',t,f(A,t));

综合实验

clear

for 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

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

综合实验一的参考放我主页里了