声明
使用R2023b版本完成,参考教材为数学实验(第三版)金正猛 王正新等编 科学出版社。参考了不少前人智慧的结晶,在此表示感谢,不保证全对,仅供大家参考,请勿抄袭,有错漏、疑问之处,烦请通过QQ或者评论区留言交流指正,欢迎大家多交流。联系QQ:3468039783
模块一:基础练习
1.1
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 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 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 我选择的是第二个实验,这是我最终版的代码,前两个版本的代码在压缩包里。 综合实验一的参考放我主页里了 精彩链接
发表评论