勤将励勉，勿望再晨。
——赠nmy
南京邮电大学通达学院《数学实验》MATLAB实验答案
- 一 声明
- 二 MATLAB下载
- 《数学实验》练习一
- 1.1
- 1.2
- 1.3
- 1.4
- 1.5
- 1.6
- 1.7
- 1.8
- 1.9
- 1.10
- 1.11
- 《数学实验》练习二
- 2.1
- 2.2
- 2.3
- 2.4
- 2.5
- 《数学实验》练习三
- 3.1
- 3.2
- 3.3
- 3.4
- 3.5
- 3.6
- 《数学实验》练习四
- 4.1
- 4.2
- 4.3
- 4.4
- 4.5
一 声明
南京邮电大学通达学院《数学实验》MATLAB实验答案
答案更新时间:2023.04.28，修改了4.2的存疑部分。已更新完成，如无错误不在更新
为了方便核算，我在代码中单独将m定义为自变量运算或者直接以m=117代入，作业中可以直接代入，即代码中不出现m。本机版本为 MATLAB R2020b
由于作者解答能力有限，难免有瑕疵错误之处，还请多多海涵！本答案仅供学习参考之用，请勿直接抄袭。有错漏之处，烦请指正。联系QQ：1415520898，如有问题可通过qq或者评论区留言方式交流。
二 MATLAB下载
这里引用@dew_142857博主的相关文章最新MATLAB R2020b超详细安装教程（附完整安装文件）实测有效，按照步骤一步步来即可，为方便同学下载，这里将文中所提向公众号索要的百度网盘链接放在下方
另外安装好的MATLAB约为96.6 GB ，请提前规划好磁盘空间。
链接：https://pan.baidu.com/s/1NExZ_v-QN4Xbu4Jk1C0dEA
提取码：7won
也可以在https://matlab.mathworks.com/注册一个账户，直接在线使用
《数学实验》练习一
1.1
log(x)——>lnx；inf——>无穷
1.2
exp(x)——>eˣ；diff(y,x,n)——>y对x的n阶导函数
1.3
第一小问答案不要忘记+C；int——>处理定积分、不定积分
1.4
2020版本
写全应该是taylor((117/200+sin(x))*cos(x),x,‘Order’,5,‘ExpansionPoint’,0)，在x=0处可省略。
2010版
1.5
本次随机的中间数据为：
[8226958330713791/9007199254740992, (2^(1/2)*469134536469018791^(1/2))/671088640, ((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2), (((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2), ((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2), (((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2), ((((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2), (((((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2), ((((((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2), (((((((((2^(1/2)*469134536469018791^(1/2))/671088640 + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2) + 117/100)^(1/2)]
1.6
本题用到的符号较多，进行下一题时使用clear清除变量
1.7
1.7.1
1.7.2
1.8
1.8.2
也可以使用下方代码，效果一样
1.9
1.10
plot是绘制二维图形，并且是x，y的表达式是已知的或者是形如y=f(x)这样确切的表达式plot函数的基本调用格式为:plot(x,y) 其中x和y为长度相同的向量，分别用于存储x坐标和y坐标数据。
ezplot是画出隐函数图形，是形如f(x,y)=0这种不能写出像y=f(x)这种函数的图形ezplot一元函数绘图函数ezplot(fun) ezplot(fun,[min,max])
fplot(y,[a,b])精确绘图
1.11
[X,Y] = meshgrid(-5:0.1:5);可以换成书上形式：
x=-5:0.1:5;y=x;
[X Y]=meshgrid(x,y);
《数学实验》练习二
2.1
第一个不动点为-0.0084
第二个不动点为119.0084
（2）先定义一个普世性的迭代方法，用M文件保存
函数收敛，只要初值不取14165^(1/2)/2+119/2 即第二个不动点，收敛值与初值的选取关系不大，总是收敛于-0.0084， 只有初值取 14165 ^(1/2)/2+119/2，迭代函数才以它为极限；
收敛值一定是不动点其一；
2.2
m=117; syms x; f=inline('1-2*abs(x-1/2)');%设定函数 x0=1/4;%设定初值 for i=1:1:10 plot(i,f(x0),'*');%用*作图，可以在括号内添加'MarkerSize',20放大点 x0=f(x0); %更新x0的值，x0类似于C语言的static类型变量 hold on %将各个点划在一张图上 end hold off
几个图像最后都是趋于0，如果没有的话要将i的终值调大，我的后三个图的i=1:1:100;
2.3
该题是P76页例二
%MARTIN函数代码 function Martin(a,b,c,N) %N为迭代次数 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))];%表示矩阵m的第n+1列。冒号表示选择所有行 end plot(m(1,:),m(2,:),'kx'); axis equal %横纵坐标采用相等单位长度 %循环迭代N次，N是预定义的数字。在循环内部，代码更新矩阵m中的值。 具体来说，该代码通过将其第一个元素设置为f(m(1，n)，m(2，n))，将其第二个元素设置为g(m(1，n))来更新m的第n列。 第一行0后面的分号表示矩阵m初始化为两行N列的列向量。m=117; Martin(m,m,m,5000) Martin(-m,-m,m,10000) Martin(-m,m/1000,-m,15000) Martin(m/1000,m/1000,0.5,20000)
2.4
（1）
%此小问无需在卷面作答，且每人选的数不一样，仔细看题目！！！ m=117; syms x; diff(subs((100*x+117)/(x^2+100),x,117^(1/3))) %对默认的变量进行一次的求导 %我取的是a=100，c=1；最后结果的绝对值应小于1才可以，否则另取 ans = 0
（2）
syms x; m=117; f=inline('(100*x+117)/(x^2+100)'); x0=10;% 任取一个初值 for i=1:20; x0=f(x0); fprintf('%g,%g\n',i,x0); end %我的运行结果 1,5.585 2,5.14893 3,4.99475 4,4.93387 5,4.90889 6,4.89849 7,4.89413 8,4.8923 9,4.89153 10,4.89121 11,4.89107 12,4.89101 13,4.89099 14,4.89098 15,4.89098 16,4.89097 17,4.89097 18,4.89097 19,4.89097 20,4.89097（3）
根据个人体会回答
函数迭代的收敛速度与初值的选取关系不大；
迭代初值对迭代的收敛性存在影响，但是这种影响存在不确定性，没有发现可循的规律；
用自己的话改一下即可
2.5
syms x; y=sin(x); y1=taylor(sin(x),x,'Order',2); y2=taylor(sin(x),x,'Order',4); y3=taylor(sin(x),x,'Order',6); fplot([y y1 y2 y3]) xlim([-3/2*pi 3/2*pi]) grid on legend('sin(x)','approximation of sin(x) up to O(x^1)','approximation of sin(x) up to O(x^3)','approximation of sin(x) up to O(x^5)')
（2）
syms x; y=sin(x); y1=taylor(sin(x),x,'Order',8); y2=taylor(sin(x),x,'Order',10); y3=taylor(sin(x),x,'Order',12); fplot([y y1 y2 y3]) xlim([-3/2*pi 3/2*pi]) grid on legend('sin(x)','approximation of sin(x) up to O(x^7)','approximation of sin(x) up to O(x^9)','approximation of sin(x) up to O(x^(11))')
(3)
《数学实验》练习三
3.1
A=str2sym('[117,117-4;6-117,10-117]');%表示符号表达式 [P,D]=eig(A); Q=inv(P); syms n; x=[1;2]; xn=P*(D.^n)*Q*x xn = (339*6^n)/2 - (337*4^n)/2 - (559*0^n)/111 2*0^n + (337*4^n)/2 - (333*6^n)/23.2
A=str2sym('[117,117-4;6-117,10-117]'); B=1/10.*A; [P,D]=eig(B); Q=inv(P); syms n; x=[1;2]; xn=P*(D.^n)*Q*x xn = (339*(3/5)^n)/2 - (337*(2/5)^n)/2 - (559*0^n)/111 2*0^n + (337*(2/5)^n)/2 - (333*(3/5)^n)/23.3
%教材P136页原题 A=[9,5;2,6]; t=[]; for i=1:20 x=2*rand(2,1)-1; t(length(t)+1,1:2)=x; for j=1:40 x=A*x; t(length(t)+1,1:2)=x; end end plot(t(:,1),t(:,2),'*') grid('on')
（2）可以看到，迭代阵列似乎在一条通过原点的直线上。
（3）
A=[9,5;2,6]; a=[]; x=2*rand(2,1)-1; for i=1:20 a(i,1:2)=x; x=A*x; end for i=1:20 if a(i,1)==0 else t=a(i,2)/a(i,1); fprintf('%g,%g\n',i,t); end end%结果 1,0.911983 2,0.551028 3,0.451391 4,0.418261 5,0.406586 6,0.402388 7,0.400867 8,0.400315 9,0.400115 10,0.400042 11,0.400015 12,0.400006 13,0.400002 14,0.400001 15,0.4 16,0.4 17,0.4 18,0.4 19,0.4 20,0.4
（4）
极限值是图像直线的斜率
按照自己语言组织下面任意一条
- 最终稳定值为迭代矩阵的特征值之一。
- 如果迭代矩阵有多个线性无关的特征向量对应于同一个特征值，那么最终稳定值将是这些特征向量线性组合的结果。
- 稳定值是迭代矩阵的特征向量，对应的特征值为1。而迭代矩阵的特征值和特征向量则可以通过特征方程来求得。
3.4
书P141相似题
m=117; A=[m-1,m;1-m,-m]; p=[0.4;0.6];%选择合适初始向量，要求和为1 [P,D]=eig(A)%P每列是特征向量，D主对角线元素是特征值 for i=1:20 p(:,i+1)=A*p(:,i); end fprintf('%2f,%2f\n',p)还可以使用下面的方法求稳定值
m=117; A=[m,1/4-m;m-3/4,1-m]; x0=[0.4;0.6]; n=10000; y = A^n * x0
结果
%A=[m,6-m;m-2,8-m] %A=[m,1/4-m;m-3/4,1-m] %A=[m-1,m;1-m,-m]
（4）ps：本题较难，可适当放弃
在线性映射迭代中，迭代矩阵的稳定性取决于其特征值的大小和分布。特征值是矩阵的一个重要性质，它描述了矩阵在线性变换下的变化情况。
如果迭代矩阵的所有特征值的绝对值都小于1，那么迭代矩阵就是稳定的，每次迭代后矩阵的元素值都会趋近于一个稳定值。
但是，如果迭代矩阵存在特征值的绝对值大于等于1，那么迭代矩阵就是不稳定的。这种情况下，每次迭代后矩阵的元素值都会趋近于无穷大或无穷小，从而导致迭代结果失效。
另外，如果迭代矩阵存在多个特征值相同的情况，那么迭代矩阵也可能不稳定。这种情况下，迭代矩阵的特征向量可能会出现非常大的幅度波动，从而导致迭代结果不可靠。
因此，对于二维矩阵的线性映射迭代，需要对迭代矩阵的特征值进行分析，以确定其稳定性。如果迭代矩阵不稳定，需要采取一些措施，如调整迭代步长或使用更稳定的迭代算法，以确保迭代结果的可靠性。
3.5
%如果默认b>a >> I=0; >> m=[]; >> n=1000; >> for a=1:n for c=a+1:n b=sqrt(c^2-a^2); if(b==floor(b))&(b>a)&(c==b+2) I=I+1;m(:,I)=[a,b,c]; end end end >> m m = 1 至 17 列 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 8 15 24 35 48 63 80 99 120 143 168 195 224 255 288 323 360 10 17 26 37 50 65 82 101 122 145 170 197 226 257 290 325 362 18 至 29 列 40 42 44 46 48 50 52 54 56 58 60 62 399 440 483 528 575 624 675 728 783 840 899 960 401 442 485 530 577 626 677 730 785 842 901 962 >> 公式：a=2m b=m^2-1 c=m^2+1(m>2,m为整数）； 即： {a,b,c}={(2u)^2,(u^2-1)^2,(u^2+1)^2}上课时默认b>a，下面给出a、b关系不确定是时的代码，无需写在试卷上
abc0=zeros(1000,3); k=0; for c=3:1000 b=c-2; a=sqrt(c^2-b^2); if(mod(a,1)==0) k=k+1; abc0(k,:)=[a b c]; end end abc=abc0(1:k,:); fprintf('所有勾股数 a b c=\n') disp(abc)3.6
for k=1:200 for b=1:999 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 1,2,8,9,18,25,32,49,50,72,81,98,121,128,162,169,200k为完全平方数或者完全平方数的二倍
预测k在[200,300]之间有200,225,242,288,289
《数学实验》练习四
4.1
% 方法一：通过法方程组求解 d0=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]; d1=sum(x);d2=sum(x.^2);b1=sum(y);b2=sum(y.*x); A=[d0,d1;d1,d2];B=[b1;b2]; u=A\B; a0=u(1) a1=u(2) error=sum((y-(a0+a1.*x)).^2) a0 = 2.8304 a1 = 4.0244 error = 0.2409%方法二：直接求解 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]; P=polyfit(x,y,1) P = 4.0244 2.8304 error=sum((y-(2.8304+4.0244.*x)).^2)%误差 error = 0.2409
4.2
（1）
%我的学号尾数是7；故数据是到1920年，对应人口是106.5，第14， %则t2=t(14)，x2=x(14) %这段代码是用来进行数据拟合的，其中变量t和x分别代表时间和数据点。代码用log函数将数据点转换成线性形式，然后使用线性回归来拟合两个数据点的斜率和截距，最后用指数函数求出x0和k，从而得到新的函数曲线。代码中的error表示新的函数曲线与原数据点的误差平方和 t=1790:10:1980; x=[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]; t1=t(1);x1=x(1); t2=t(14);x2=x(14); %此步根据学号不同而不同 A=[1,t1;1,t2]; b=[log(x1);log(x2)]; u=A\b; x0=exp(u(1)) k=u(2) error=sum((x0*exp(k*t)-x).^2) x0 = 6.5242e-20 k = 0.0254 error = 1.2278e+05（2）
t=1790:10:1920; x=[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]; %我的数据是到1920，所以上面的数据是截到1920年对应的106.5 y=log(x); m=length(t); A=[m,sum(t);sum(t),sum(t.^2)]; b=[sum(y);y*t']; u=A\b; x0=exp(u(1)) k=u(2) error=sum((x0*exp(k*t)-x).^2) x0 = 2.7207e-20 k = 0.0260 error = 681.95884.3
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.01]; f=@(a,x)a(1)./(1+a(2)*exp(-a(3)*x)); [A,resnorm]=lsqcurvefit(f,a,x,y) f(A,20) A = 1.0e+03 * 5.7882 0.0025 0.0001 resnorm = 3.3995e+04 ans = 4.7438e+034.4
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=27; while f(A,t+1)-f(A,t)>=10 t=t+1; end f(A,t) A = 1.0e+03 * 5.3860 0.0021 0.0001 0.0000 resnorm = 9.1025e+03 ans = 5.3409e+034.5
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=[2,0.1,0.1];%r、k、a f=@(a,x)a(1)*exp(a(2)*x+a(3)); [A,resnorm]=lsqcurvefit(f,a,x,y) t=27; while f(A,t+1)-f(A,t)>=10 t=t+1; end f(A,t) A = 2.4511 0.3152 -0.0841 resnorm = 2.3628e+08 ans = Inf
![[工业自动化-1]：PLC架构与工作原理](https://img-blog.csdnimg.cn/ce10a1471ed14382bc58364cf8bd5209.png)
还没有评论，来说两句吧...