MATLAB习题答案(清华大学)

高等应用数学问题MATLAB求解习题参考解答（薛定宇著）目录第1章计算机数学语言概述2第2章MATLAB语言程序设计基础5第3章微积分问题的计算机求解17第4章线性代数问题的计算机求解29第5章积分变换与复变函数问题的计算机求解43第6章代数方程与最优化问题的计算机求解53第7章微分方程问题的计算机求解71第8章数据插值、函数逼近问题的计算机求解93第9章概率论与数理统计问题的计算机求解114第10章数学问题的非传统解法127第A章自由数学语言Scilab简介136第1章计算机数学语言概述1在你的机器上安装MATLAB语言环境，并键入demo命令，由给出的菜单系统和对话框原型演示程序，领略MATLAB语言在求解数学问题方面的能力与方法。【求解】在MATLAB提示符>>下键入demo命令，则将打开如图1-1所示的窗口，窗口左侧的列表框可以选择各种不同组合的演示内容。图1-1MATLAB演示程序界面例如，用户选择MATLAB!Graphics!VolumeVlsulization演示，则将得出如图1-2所示的演示说明，单击其中的Runthisdemo栏目，则将得出如图1-3所示的演示界面。用户可以在该界面下按按钮，逐步演示相关内容，而实现这样演示的语句将在该程序界面的下部窗口中给出。2作者用MATLAB语言编写了给出例子的源程序，读者可以自己用type语句阅读一下源程序，对照数学问题初步理解语句的含义，编写的源程序说明由下表列出。第1章计算机数学语言概述3图1-2MATLAB演示程序界面举例序号文件名程序说明例1.1c1ex1.m利用MATLAB的符号运算工具箱求解微分问题例1.2c1ex2.m分别利用MATLAB的符号运算工具箱和数值运算功能求解多项式方程，其中用数值方法得出的结果有误差例1.3c1ex3.m分别利用MATLAB的符号运算工具箱和数值运算功能计算Hilbert矩阵的行列式，其中用数值方法得出的结果有很大误差例1.4c1ex4.m令x1=y;x2=y_，则可以将原来的二阶微分方程转换成一阶微分方程组，然后就可以求解微分方程的数值解了，原方程是非线性微分方程，故不存在解析解。ode45()函数可以求解常微分方程组，而dde23()可以求解延迟微分方程，或更直观地采用Simulink绘制求解框图。例1.5c1ex5.m线性规划问题调用最优化工具箱中的linprog()函数可以立即得出结果，若想求解整数规划问题，则需要首先安装整数规划程序ipslvmex()。4第1章计算机数学语言概述图1-3MATLAB体视化演示程序界面第2章MATLAB语言程序设计基础1启动MATLAB环境，并给出语句tic,A=rand(500);B=inv(A);norm(A*B-eye(500)),toc，试运行该语句，观察得出的结果，并利用help命令对你不熟悉的语句进行帮助信息查询，逐条给出上述程序段与结果的解释。【求解】在MATLAB环境中感触如下语句，则可以看出，求解500£500随机矩阵的逆，并求出得出的逆矩阵与原矩阵的乘积，得出和单位矩阵的差，得出范数。一般来说，这样得出的逆矩阵精度可以达到10¡12。>>tic,A=rand(500);B=inv(A);norm(A*B-eye(500)),tocans=1.2333e-012Elapsedtimeis1.301000seconds.2试用符号元素工具箱支持的方式表达多项式f(x)=x5+3x4+4x3+2x2+3x+6，并令x=s¡1s+1，将f(x)替换成s的函数。【求解】可以先定义出f函数，则由subs()函数将x替换成s的函数>>symssxf=x^5+3*x^4+4*x^3+2*x^2+3*x+6;F=subs(f,x,(s-1)/(s+1))F=(s-1)^5/(s+1)^5+3*(s-1)^4/(s+1)^4+4*(s-1)^3/(s+1)^3+2*(s-1)^2/(s+1)^2+3*(s-1)/(s+1)+63用MATLAB语句输入矩阵A和B矩阵①A=266412344321234132413775;②B=26641+4j2+3j3+2j4+1j4+1j3+2j2+3j1+4j2+3j3+2j4+1j1+4j3+2j2+3j4+1j1+4j3775前面给出的是4£4矩阵，如果给出A(5;6)=5命令将得出什么结果？【求解】用课程介绍的方法可以直接输入这两个矩阵>>A=[1234;4321;2341;3241]A=12346第2章MATLAB语言程序设计基础432123413241若给出A(5,6)=5命令，虽然这时的行和列数均大于B矩阵当前的维数，但仍然可以执行该语句，得出>>A(5,6)=5A=123400432100234100324100000005复数矩阵也可以用直观的语句输入>>B=[1+4i2+3i3+2i4+1i;4+1i3+2i2+3i1+4i;2+3i3+2i4+1i1+4i;3+2i2+3i4+1i1+4i];B=1.0000+4.0000i2.0000+3.0000i3.0000+2.0000i4.0000+1.0000i4.0000+1.0000i3.0000+2.0000i2.0000+3.0000i1.0000+4.0000i2.0000+3.0000i3.0000+2.0000i4.0000+1.0000i1.0000+4.0000i3.0000+2.0000i2.0000+3.0000i4.0000+1.0000i1.0000+4.0000i4假设已知矩阵A，试给出相应的MATLAB命令，将其全部偶数行提取出来，赋给B矩阵，用A=magic(8)命令生成A矩阵，用上述的命令检验一下结果是不是正确。【求解】魔方矩阵可以采用magic()生成，子矩阵也可以提取出来>>A=magic(8),B=A(2:2:end,:)A=64236160675795554121351501617474620214342244026273736303133323435292838392541232244451918484915145253111056858595462631B=955541213515016第2章MATLAB语言程序设计基础7402627373630313341232244451918488585954626315用MATLAB语言实现下面的分段函数y=f(x)=8<:h;x>Dh=Dx;jxj6D¡h;x<¡D。【求解】两种方法，其一，巧用比较表达式解决>>y=h*(x>D)+h/D*x.*(abs(x)<=D)-h*(x<-D);另外一种方法，用循环语句和条件转移语句>>fori=1:length(x)ifx(i)>D,y(i)=h;elseifabs(x(i))<=D,y(i)=h/D*x(i);else,y(i)=-h;endend其中，前者语句结构简单，但适用范围更广，允许使用矩阵型x，后者只能使用向量型的x，但不能处理矩阵问题。6用数值方法可以求出S=X63i=02i=1+2+4+8+¢¢¢+262+263，试不采用循环的形式求出和式的数值解。由于数值方法采用double形式进行计算的，难以保证有效位数字，所以结果不一定精确。试采用符号运算的方法求该和式的精确值。【求解】用符号运算的方式可以采用下面语句>>sum(sym(2).^[1:63])ans=18446744073709551614由于结果有19位数值，所以用double型不能精确表示结果，该数据类型最多表示16位有效数字。其实用符号运算方式可以任意保留有效数字，例如可以求200项的和或1000项的和可以由下面语句立即得出。>>sum(sym(2).^[1:200])ans=3213876088517980551083924184682325205044405987565585670602750>>sum(sym(2).^[1:1000])ans=2143017214372534641896850098120003621122809623411067214887500776740702102249872244986396757631391716255189345835106293650374290571384628088第2章MATLAB语言程序设计基础719691551493971496078691355496484619708421492101247422837559083643060929499671638825347975351183310878921541258291423929553730843353208596633052487736744113361387507编写一个矩阵相加函数matadd()，使其具体的调用格式为A=matadd(A1,A2,A3,¢¢¢)，要求该函数能接受任意多个矩阵进行解法运算。【求解】可以编写下面的函数，用varargin变量来表示可变输入变量functionA=mat_add(varargin)A=0;fori=1:length(varargin),A=A+varargin{i};end如果想得到合适的错误显示，则可以试用try,catch结构。functionA=mat_add(varargin)tryA=0;fori=1:length(varargin),A=A+varargin{i};endcatch,error(lasterr);end8自己编写一个MATLAB函数，使它能自动生成一个m£m的Hankel矩阵，并使其调用格式为v=[h1;h2;hm;hm+1;¢¢¢;h2m¡1];H=myhankel(v)【求解】解决这样的问题可以有多种方法：①最直接的方法，Hi;j=hi+j¡1，利用双重循环functionH=myhankel(v)m=(length(v)+1)/2;%严格说来还应该判定给定输入向量长度奇偶性fori=1:m,forj=1:mH(i,j)=v(i+j-1);end,end②考虑某一行(或列)，ai=[hi;hi+1;¢¢¢;hi+m¡1]，就可以用单重循环生成Hankel矩阵了functionH=myhankel(v)m=(length(v)+1)/2;%严格说来还应该判定给定输入向量长度奇偶性fori=1:m,H(i,:)=v(i:i+m-1);end③利用现有的hankel()函数，则functionH=myhankel(v)m=(length(v)+1)/2;%严格说来还应该判定给定输入向量长度奇偶性H=hankel(v(1:m),v(m:end));第2章MATLAB语言程序设计基础99已知Fibonacci数列由式ak=ak¡1+ak¡2;k=3;4;¢¢¢可以生成，其中初值为a1=a2=1，试编写出生成某项Fibonacci数值的MATLAB函数，要求①函数格式为y=fib(k)，给出k即能求出第k项ak并赋给y向量；②编写适当语句，对输入输出变量进行检验，确保函数能正确调用；③利用递归调用的方式编写此函数。【求解】假设fib(n)可以求出Fibonacci数列的第n项，所以对n>3则可以用k=fib(n¡1)+fib(n¡2)可以求出数列的n+1项，这可以使用递归调用的功能，而递归调用的出口为1。综上，可以编写出M-函数。functiony=fib(n)ifround(n)==n&n>=1ifn>=3y=fib(n-1)+fib(n-2);else,y=1;endelseerror('nmustbepositiveinteger.')end例如，n=10可以求出相应的项为>>fib(10)ans=55现在需要比较一下递归实现的速度和循环实现的速度>>tic,fib(20),tocans=832040elapsed_time=62.0490>>tic,a=[11];fori=3:30,a(i)=a(i-1)+a(i-2);end,a(30),tocans=832040elapsed_time=0.0100应该指出，递归的调用方式速度较慢，比循环语句慢很多，所以不是特别需要，解这样问题没有必要用递归调用的方式。10第2章MATLAB语言程序设计基础10由矩阵理论可知，如果一个矩阵M可以写成M=A+BCBT,并且其中A,B,C为相应阶数的矩阵，则M矩阵的逆矩阵可以由下面的算法求出M¡1=³A+BCBT´¡1=A¡1¡A¡1B³C¡1+BTA¡1B´¡1BTA¡1试根据上面的算法用MATLAB语句编写一个函数对矩阵M进行求逆，并通过一个小例子来检验该程序，并和直接求逆方法进行精度上的比较。【求解】编写这个函数functionMinv=part_inv(A,B,C)Minv=inv(A)-inv(A)*B*inv(inv(C)+B'*inv(A)*B)*B'*inv(A);假设矩阵为M=2664515036165077603236608748163248683775且已知该矩阵可以分解成A=266410000200003000043775;B=266412342340340040003775;C=266440000300002000013775对这个例子。可以>>M=[51503616;50776032;36608748;16324868];iM=inv(M);%数值逆，直接解法iM=0.0553-0.03890.00170.0041-0.03890.0555-0.0210-0.00210.0017-0.02100.0328-0.01370.0041-0.0021-0.01370.0244>>A=diag([1234]);B=hankel([1234]);C=diag([4321]);iM1=part_inv(A,B,C)%分块矩阵的求解方法iM1=0.0553-0.03890.00170.0041-0.03890.0555-0.0210-0.00210.0017-0.02100.0328-0.01370.0041-0.0021-0.01370.0244乍看结果，似乎二者完全一致，实际上数值算法是有区别的。我们这里用解析方法得出矩阵的逆，然后用下面的语句比较两个结果的精度第2章MATLAB语言程序设计基础11>>M1=sym(M);iM0=inv(M1)iM0=[10713/193751,-7546/193751,332/193751,796/193751][-7546/193751,10759/193751,-4068/193751,-416/193751][332/193751,-4068/193751,19075/581253,-2652/193751][796/193751,-416/193751,-2652/193751,18919/775004]>>norm(double(iM0)-iM)%直接求解的误差范数ans=2.7990e-017>>norm(double(iM0)-iM1)%间接求解的误差范数ans=3.6583e-016可见，用间接方法得出的逆矩阵误差更大，因为在这里新编写的函数中inv()函数使用了多次，由此产生很大的传递误差。由此可以得出结论：如果某问题存在直接解，则尽量别使用间接方法，以加大传递误差。11下面给出了一个迭代模型(xk+1=1+yk¡1:4x2kyk+1=0:3xk写出求解该模型的M-函数，如果取迭代初值为x0=0,y0=0，那么请进行30000次迭代求出一组x和y向量，然后在所有的xk和yk坐标处点亮一个点(注意不要连线)，最后绘制出所需的图形。提示这样绘制出的图形又称为Henon引力线图，它将迭代出来的随机点吸引到一起，最后得出貌似连贯的引力线图。【求解】用循环形式解决此问题，可以得出如图2-1所示的Henon引力线图。>>x=0;y=0;fori=1:29999x(i+1)=1+y(i)-1.4*x(i)^2;y(i+1)=0.3*x(i);endplot(x,y,'.')上述的算法由于动态定义x和y，所以每循环一步需要重新定维，这样做是很消耗时间的，所以为加快速度，可以考虑预先定义这两个变量，如给出x=zeros(1,30000)。12用MATLAB语言的基本语句显然可以立即绘制一个正三角形，试结合循环结构，编写一个小程序，在同一个坐标系下绘制出该正三角形绕其中心旋转后得出的一系列三角形，还可以调整旋转步距观察效果。12第2章MATLAB语言程序设计基础−1.5−1−0.500.511.5−0.4−0.3−0.2−0.100.10.20.30.4图2-1Henon引力线图【求解】假设正三角形逆时针旋转µ度，则可以得出如图2-2a所示的示意图，三角形的三个顶点为(cosµ;sinµ),(cos(µ+120±);sin(µ+120±)),(cos(µ+240±);sin(µ+240±))，可以绘制出其曲线，如图2-2b所示，试减小步距，如选择¢µ=2;1;0:1，观察效果。.............................................................µxy-6(a)示意图−1−0.500.51−1−0.500.51(b)曲线绘制效果图2-2曲线绘制>>t=[0,120,240,0]*pi/180;%变换成弧度xxx=[];yyy=[];fori=0:5:360tt=i*pi/180;xxx=[xxx;cos(tt+t)];yyy=[yyy;sin(tt+t)];end第2章MATLAB语言程序设计基础13plot(xxx',yyy','r'),axis('square')13选择合适的步距绘制出下面的图形sinµ1t¶，其中t2(¡1;1)。【求解】用普通的绘图形式，选择等间距，得出如图2-3a所示的曲线，其中x=0左右显得粗糙。>>t=-1:0.03:1;y=sin(1./t);plot(t,y)选择不等间距方法，可以得出如图2-3b所示的曲线。>>t=[-1:0.03:-0.25,-0.248:0.001:0.248,0.25:.03:1];y=sin(1./t);plot(t,y)−1−0.500.51−1−0.500.51(a)等间距曲线绘制−1−0.500.51−1−0.500.51(b)不等间距曲线绘制图2-3不同自变量选取下的sin(1=t)曲线14对合适的µ范围选取分别绘制出下列极坐标图形①½=1:0013µ2，②½=cos(7µ=2)，③½=sin(µ)=µ，④½=1¡cos3(7µ)【求解】绘制极坐标曲线的方法很简单，用polar(µ,½)即可以绘制出极坐标图，如图2-4所示。注意绘制图形时的点运算：>>t=0:0.01:2*pi;subplot(221),polar(t,1.0013*t.^2),%(a)subplot(222),t1=0:0.01:4*pi;polar(t1,cos(7*t1/2))%(b)subplot(223),polar(t,sin(t)./t)%(c)subplot(224),polar(t,1-(cos(7*t)).^3)15用图解的方式找到下面两个方程构成的联立方程的近似解。x2+y2=3xy2;x3¡x2=y2¡y【求解】这两个方程应该用隐式方程绘制函数ezplot()来绘制，交点即方程的解，如图2-5a所示。14第2章MATLAB语言程序设计基础204030210602409027012030015033018000.5130210602409027012030015033018000.513021060240902701203001503301800123021060240902701203001503301800图2-4极坐标图>>ezplot('x^2+y^2-3*x*y^2');holdonezplot('x^3-x^2=y^2-y')可用局部放大的方法求出更精确的值，如图2-5b所示。从图上可以精确读出两个交点，(0:4012;¡0:8916)，(1:5894;0:8185)。试将这两个点分别代入原始方程进行验证。−6−4−20246−6−4−20246xyx3−x2=y2−y=0(a)两个方程的曲线，交点为解1.58941.58941.58941.58941.58941.58941.58940.81850.81850.81850.81850.8185xyx3−x2=y2−y=0(b)局部放大区域图2-5二元联立方程的图解法第2章MATLAB语言程序设计基础1516请分别绘制出xy和sin(xy)的三维图和等高线。【求解】(a)给出下面命令即可，得出的图形如图2-6a、b所示。>>[x,y]=meshgrid(-1:.1:1);surf(x,y,x.*y),figure;contour(x,y,x.*y,30)(b)给出下面命令即可，得出的图形如图2-6c、d所示。>>[x,y]=meshgrid(-pi:.1:pi);surf(x,y,sin(x.*y)),figure;contour(x,y,sin(x.*y),30)−101−101−1−0.500.51(a)xy三维图−1−0.500.51−1−0.500.51(b)xy等高线−505−505−1−0.500.51(c)sin(xy)三维图−3−2−10123−3−2−10123(d)sin(xy)等高线图2-6三维图与等高线17在图形绘制语句中，若函数值为不定式NaN，则相应的部分不绘制出来，试利用该规律绘制z=sinxy的表面图，并剪切下x2+y260:52的部分。【求解】给出下面命令可以得出矩形区域的函数值，再找出x2+y260:52区域的坐标，将其函数值设置成NaN，最终得出如图2-7所示的曲面。16第2章MATLAB语言程序设计基础>>[x,y]=meshgrid(-1:.1:1);z=sin(x.*y);ii=find(x.^2+y.^2<=0.5^2);z(ii)=NaN;surf(x,y,z)−1−0.500.51−1−0.500.51−1−0.500.51图2-7得出的三维图第3章微积分问题的计算机求解1试求出如下极限。①limx!1(3x+9x)1x，②limx!1(x+2)x+2(x+3)x+3(x+5)2x+5【求解】极限问题由下面的语句可以直接求出。>>symsx;f=(3^x+9^x)^(1/x);limit(f,x,inf)ans=9>>symsx;f=(x+2)^(x+2)*(x+3)^(x+3)/(x+5)^(2*x+5);limit(f,x,inf)ans=exp(-5)2试求下面的双重极限。①limx!¡1y!2x2y+xy3(x+y)3，②limx!0y!0xypxy+1¡1，③limx!0y!01¡cos¡x2+y2¢¡x2+y2¢ex2+y2。【求解】双重极限问题可以由下面语句直接求解。>>symsxyfa=(x^2*y+x*y^3)/(x+y)^3;limit(limit(fa,x,-1),y,2)ans=-6>>fb=x*y/(sqrt(x*y+1)-1);limit(limit(fb,x,0),y,0)ans=2>>fc=(1-cos(x^2+y^2))*exp(x^2+y^2)/(x^2+y^2);limit(limit(fc,x,0),y,0)ans=03求出下面函数的导数。①y(x)=qxsinxp1¡ex，②y=1¡pcosaxx(1¡cospax)③atanyx=ln(x2+y2)，④y(x)=¡1nalnxn+axn;n>018第3章微积分问题的计算机求解【求解】由求导函数diff()可以直接得出如下结果，其中③为隐函数，故需要用隐函数求导公式得出导数。>>symsx;f=sqrt(x*sin(x)*sqrt(1-exp(x)));simple(diff(f))ans=1/2/(x*sin(x)*(1-exp(x))^(1/2))^(1/2)*(sin(x)*(1-exp(x))^(1/2)+x*cos(x)*(1-exp(x))^(1/2)-1/2*x*sin(x)/(1-exp(x))^(1/2)*exp(x))>>symsaxy=(1-sqrt(cos(a*x)))/(x*(1-cos(sqrt(a*x))))simple(diff(y))ans=1/2/cos(a*x)^(1/2)*sin(a*x)*a/x/(1-cos((a*x)^(1/2)))-(1-cos(a*x)^(1/2))/x^2/(1-cos((a*x)^(1/2)))-1/2*(1-cos(a*x)^(1/2))/x/(1-cos((a*x)^(1/2)))^2*sin((a*x)^(1/2))/(a*x)^(1/2)*a>>f=atan(y/x)-log(x^2+y^2);f1=simple(-diff(f,x)/diff(f,y))f1=(y+2*x)/(x-2*y)>>symsnpositive;symsa;f=-log((x^n+a)/x^n)/(n*a);diff(f,x)ans=-(n/x-(x^n+a)/(x^n)*n/x)/(x^n+a)*x^n/n/a用LATEX表示上面的结果，则①1=2µsin(x)p1¡ex+xcos(x)p1¡ex¡1=2xsin(x)exp1¡ex¶1qxsin(x)p1¡ex②1=2sin(ax)apcos(ax)x(1¡cos(pax))¡1¡pcos(ax)x2(1¡cos(pax))¡1=2³1¡pcos(ax)´sin(pax)ax(1¡cos(pax))2pax③y+2xx¡2y④¡µnx¡(xn+a)nxnx¶xn(xn+a)¡1n¡1a¡4试求出y(t)=s(x¡1)(x¡2)(x¡3)(x¡4)函数的4阶导数。【求解】高阶导数可以由下面语句直接得出>>symsaxf=sqrt((x-1)*(x-2)/(x-3)/(x-4));simple(diff(f,x,4))ans=第3章微积分问题的计算机求解193*(16*x^11-392*x^10+4312*x^9-28140*x^8+121344*x^7-364560*x^6+783552*x^5-1214604*x^4+1342560*x^3-1015348*x^2+474596*x-103741)/((x-1)*(x-2)/(x-3)/(x-4))^(7/2)/(x-3)^8/(x-4)^83µ16x11¡392x10+4312x9¡28140x8+121344x7¡364560x6+783552x5¡1214604x4+1342560x3¡1015348x2+474596x¡103741¶µ(x¡1)(x¡2)(x¡3)(x¡4)¶7=2(x¡3)8(x¡4)85在高等数学中，求解分子和分母均同时为0或1时，分式极限时可使用L'H^opital法则，即对分子分母分别求导数，再由比值得出，试用该法则limx!0ln(1+x)ln(1¡x)¡ln(1¡x2)x4，并和直接求出的极限结果相比较。【求解】从给出的分母看，若想使之在x=0处的值不为0，则应该对其求4阶导数，同样，还应该对分子求4阶导数，将x=0代入结果，这样就可以使用L'H^opital法则求出极限了。>>symsx;n=log(1+x)*log(1-x)-log(1-x^2);d=x^4;n4=diff(n,x,4);d4=diff(d,x,4);n4=subs(n4,x,0);L=n4/d4L=1/12现在直接求极限可以验证上述结果是正确的。>>limit(n/d,x,0)ans=1/126已知参数方程½x=lncosty=cost¡tsint，试求出dydx和d2ydx2¯¯¯¯t=p=3。【求解】参数方程的导数可以由下面语句直接求出。>>symst;x=log(cos(t));y=cos(t)-t*sin(t);diff(y,t)/diff(x,t)ans=-(-2*sin(t)-t*cos(t))/sin(t)*cos(t)>>f=diff(y,t,2)/diff(x,t,2);subs(f,t,sym(pi)/3)ans=3/8-1/24*pi*3^(1/2)7假设u=cos¡1rxy，试验证@2u@x@y=@2u@y@x。【求解】证明二者相等亦可以由二者之差为零来证明，故由下面的语句直接证明。20第3章微积分问题的计算机求解>>symsxy;u=acos(x/y);diff(diff(u,x),y)-diff(diff(u,y),x)ans=08设(xu+yv=0yu+xv=1，试求解@2u@x@y。【求解】用下面的语句可以直接得出如下结果。>>symsxyuv[u,v]=solve('x*u+y*v=0','y*u+x*v=1','u,v');diff(diff(u,x),y)ans=2/(x^2-y^2)^2*x+8*y^2/(x^2-y^2)^3*x9假设f(x;y)=Zxy0e¡t2dt，试求xy@2@x2¡2@2@x@y+@2@y2。【求解】由下面的命令可以得出所需结果。>>symsxytf=int(exp(-t^2),t,0,x*y);x/y*diff(f,x,2)-2*diff(diff(f,x),y)+diff(f,y,2)simple(ans)ans=-2*exp(-x^2*y^2)*(-x^2*y^2+1+x^3*y)10假设已知函数矩阵f(x;y;z)=·3x+eyzx3+y2sinz¸，试求出其Jacobi矩阵。【求解】Jacobi矩阵可以由下面的语句直接得出。>>symsxyzF=[3*x+exp(y)*z;x^3+y^2*sin(z)];jacobian(F,[x,y,z])ans=[3,exp(y)*z,exp(y)][3*x^2,2*y*sin(z),y^2*cos(z)]11试求解下面的不定积分问题。①I(x)=¡Z3x2+ax2(x2+a)2dx，②I(x)=Zpx(x+1)px+p1+xdx③I(x)=Zxeaxcosbxdx，④I(t)=Zeaxsinbxsincxdx【求解】①该不定积分可以由下面的命令直接求出第3章微积分问题的计算机求解21>>symsxaf=(3*x^2+a)/(x^2+(x^2+a)^2);int(f,x)ans=12/(4+16*a)/(2+4*a+2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a+2*(1+4*a)^(1/2))^(1/2))+48/(4+16*a)/(2+4*a+2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a+2*(1+4*a)^(1/2))^(1/2))*a+12/(4+16*a)/(2+4*a+2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a+2*(1+4*a)^(1/2))^(1/2))*(1+4*a)^(1/2)+16/(4+16*a)/(2+4*a+2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a+2*(1+4*a)^(1/2))^(1/2))*(1+4*a)^(1/2)*a+12/(4+16*a)/(2+4*a-2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a-2*(1+4*a)^(1/2))^(1/2))+48/(4+16*a)/(2+4*a-2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a-2*(1+4*a)^(1/2))^(1/2))*a-12/(4+16*a)/(2+4*a-2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a-2*(1+4*a)^(1/2))^(1/2))*(1+4*a)^(1/2)-16/(4+16*a)/(2+4*a-2*(1+4*a)^(1/2))^(1/2)*atan(2*x/(2+4*a-2*(1+4*a)^(1/2))^(1/2))*(1+4*a)^(1/2)*a②可以用下面的语句求出问题的解>>symsx;f=sqrt(x*(x+1))/(sqrt(x)+sqrt(x+1));int(f,x);latex(ans)并将其显示如下2=15px(x+1)x(3x+5)px+1¡2=15px(x+1)(x+1)(¡2+3x)px③可以求出下面的结果>>symsabxf=x*exp(a*x)*cos(b*x);int(f,x);latex(ans)其数学显示为Ãaxa2+b2¡a2¡b2(a2+b2)2!eaxcos(bx)¡Ã¡bxa2+b2+2a(a2+b2)2!eaxsin(bx)④用下面的语句求解，得>>symsxabc;f=exp(a*x)*sin(b*x)*sin(c*x);latex(int(f,x))22第3章微积分问题的计算机求解亦即1=2aeaxcos((b¡c)xa2+(b¡c)2¡1=2(¡b+c)eaxsin((b¡c)x)a2+(b¡c)2¡1=2aeaxcos((b+c)xa2+(b+c)2+1=2(¡b¡c)eaxsin((b+c)x)a2+(b+c)212试求出下面的定积分或无穷积分。①I=Z10cosxpxdx，②I=Z101+x21+x4dx【求解】①可以直接求解>>symsx;int(cos(x)/sqrt(x),x,0,inf)ans=1/2*2^(1/2)*pi^(1/2)②可以得出>>symsx;int((1+x^2)/(1+x^4),x,0,1)ans=1/4*2^(1/2)*pi13假设f(x)=e¡5xsin(3x+p=3)，试求出积分函数R(t)=Zt0f(x)f(t+x)dx。【求解】定义了x的函数，则可以由subs()函数定义出t+x的函数，这样由下面的语句可以直接得出R函数。>>symsxt;f=exp(-5*x)*sin(3*x+sym(pi)/3);R=int(f*subs(f,x,t+x),x,0,t);simple(R)ans=1/1360*(15*exp(t)^10*3^(1/2)*cos(3*t)-25*cos(9*t)+25*exp(t)^10*3^(1/2)*sin(3*t)-68*cos(3*t)-15*3^(1/2)*cos(9*t)-25*3^(1/2)*sin(9*t)-15*exp(t)^10*sin(3*t)+15*sin(9*t)+93*exp(t)^10*cos(3*t))/exp(t)^15该结果可以写成1136015(et)10p3cos(3t)¡68cos(3t)¡15(et)10sin(3t)¡25p3sin(9t)+25(et)10p3sin(3t)+15sin(9t)¡25cos(9t)¡15p3cos(9t)+93(et)10cos(3t)(et)1514对a的不同取值试求出I=Z10cosax1+x2dx。【求解】由下面的循环结构可以得出不同a值下的无穷积分值，并可以绘制出它们之间关系的曲线，如图3-1所示。第3章微积分问题的计算机求解23>>symsxa;f=cos(a*x)/(1+x^2);aa=[0:0.1:pi];y=[];forn=aab=int(subs(f,a,n),x,0,inf);y=[y,double(b)];endplot(aa,y)00.511.522.533.500.20.40.60.811.21.41.6图3-1不同a值下的积分值曲线15试对下面函数进行Fourier幂级数展开。①f(x)=(p¡jxj)sinx;¡p6x<p②f(x)=ejxj;¡p6x<p③f(x)=(2x=l;0<x<l=22(l¡x)=l;l=2<x<l，且l=p。【求解】①可以立即由下面的语句求出。>>symsx;f=(sym(pi)-abs(x))*sin(x);[A,B,F]=fseries(f,x,10,-pi,pi);FF=1/2*pi*sin(x)+16/9/pi*sin(2*x)+32/225/pi*sin(4*x)+48/1225/pi*sin(6*x)+64/3969/pi*sin(8*x)+80/9801/pi*sin(10*x)该结果在LATEX下可以显示为12psinx+169sin2xp+32225sin4xp+481225sin6xp+643969sin8xp+809801sin10xp②可以由下面语句求解，并得出数学公式为>>symsx;f=exp(abs(x));[A,B,F]=fseries(f,x,10,-pi,pi);F得出的解析解为1=22ep¡2p+(¡ep¡1)cos(x)p+(2=5ep¡2=5)cos(2x)p+(¡1=5ep¡1=5)cos(3x)p24第3章微积分问题的计算机求解+(2=17ep¡2=17)cos(4x)p+(¡1=13ep¡1=13)cos(5x)p+¡237ep¡237¢cos(6x)p+(¡1=25ep¡1=25)cos(7x)p+¡265ep¡265¢cos(8x)p+(¡1=41ep¡1=41)cos(9x)p+¡2101ep¡2101¢cos(10x)p进一步观察结果可见，该式子可以手工化简，例如提取系数(ep¡(¡1)n)=p。或对各项系数逐项求值(保留10位有效数字)>>vpa(F,10)ans=7.047601355-7.684221126*cos(x)+2.819040541*cos(2.*x)-1.536844225*cos(3.*x)+.8291295709*cos(4.*x)-.5910939328*cos(5.*x)+.3809514246*cos(6.*x)-.3073688450*cos(7.*x)+.2168492724*cos(8.*x)-.1874200274*cos(9.*x)+.1395564625*cos(10.*x)③似乎求解起来更困难，巧妙利用符号运算工具箱中的heaviside()函数，则可以将原函数表示成f(x)=2¤heaviside³x¡p2´¡2pxjx¡p=2jx¡p=2这样就可以用下面的语句求出函数的Fourier级数。>>symsx;pi1=sym(pi);f=2*heaviside(x-pi1/2)-2/pi1*x*abs(x-pi1/2)/(x-pi1/2);[a,b,F]=fseries(f,x,10,-pi,pi);FF=-1/4+4/pi^2*cos(x)+(4/pi+2)/pi*sin(x)-2/pi^2*cos(2*x)-1/pi*sin(2*x)+4/9/pi^2*cos(3*x)+(-4/9/pi+2/3)/pi*sin(3*x)-1/2/pi*sin(4*x)+4/25/pi^2*cos(5*x)+(4/25/pi+2/5)/pi*sin(5*x)-2/9/pi^2*cos(6*x)-1/3/pi*sin(6*x)+4/49/pi^2*cos(7*x)+(-4/49/pi+2/7)/pi*sin(7*x)-1/4/pi*sin(8*x)+4/81/pi^2*cos(9*x)+(4/81/pi+2/9)/pi*sin(9*x)-2/25/pi^2*cos(10*x)-1/5/pi*sin(10*x)¡1=4+4cos(x)p2+¡4p¡1+2¢sin(x)p¡2cos(2x)p2¡sin(2x)p+4=9cos(3x)p2+¡¡4=9p¡1+2=3¢sin(3x)p¡1=2sin(4x)p+425cos(5x)p2+¡425p¡1+2=5¢sin(5x)p¡2=9cos(6x)p2¡1=3sin(6x)p+449cos(7x)p2+¡¡449p¡1+2=7¢sin(7x)p¡1=4sin(8x)p+481cos(9x)p2+¡481p¡1+2=9¢sin(9x)p¡225cos(10x)p2¡1=5sin(10x)p16试求出下面函数的Taylor幂级数展开。第3章微积分问题的计算机求解25①Zx0sinttdt②lnµ1+x1¡x¶③ln³x+p1+x2´④(1+4:2x2)0:2⑤e¡5xsin(3x+p=3)分别关于x=0、x=a的幂级数展开⑥对f(x;y)=1¡cos¡x2+y2¢¡x2+y2¢ex2+y2关于x=1;y=0进行二维Taylor幂级数展开。【求解】由下面的语句可以分别求出各个函数的幂级数展开，由latex(ans)函数可以得出下面的数学表示形式。>>symstx;f=int(sin(t)/t,t,0,x);taylor(f,x,15)>>symsx;f=log((1+x)/(1-x)),taylor(f,x,15)>>symsx;f=log(x+sqrt(1+x^2));taylor(f,x,15)>>symsx;f=(1+4.2*x^2)^0.2;taylor(f,x,13)①x¡1=18x3+1600x5¡135280x7+13265920x9¡1439084800x11+180951270400x13②2x+2=3x3+2=5x5+2=7x7+2=9x9+2=11x11+2=13x13③x¡1=6x3+340x5¡5112x7+351152x9¡632816x11+23113312x13④1+2125x2¡882625x4+5556615625x6¡4084101390625x8+162955629948828125x10¡1368827291161220703125x12⑤该函数的前4项展开为>>symsxa;f=exp(-5*x)*sin(3*x+sym(pi)/3);taylor(f,x,4,a)e¡5a³3a+3´+³3e¡5a³3a+3´¡5e¡5a³3a+3´´(x¡a)+³8e¡5a³3a+3´¡15e¡5a³3a+3´´(x¡a)2+³33e¡5a³3a+3´+5=3e¡5a³3a+3´´(x¡a)3⑥该函数需要使用Maple的展开函数。>>symsxy;f=(1-cos(x^2+y^2))/((x^2+y^2)*exp(x^2+y^2));F=maple('mtaylor',f,'[x=1,y]',4)1¡cos(1)e1+(2sin(1)¡4+4cos(1))(x¡1)e1+(¡6cos(1)¡7sin(1)+8)(x¡1)2e1+(sin(1)¡2+2cos(1))y2e1+¡¡343+323sin(1)+16=3cos(1)¢(x¡1)3e1+(¡8cos(1)+10¡8sin(1))y2(x¡1)e1+¡836¡6cos(1)¡343sin(1)¢(x¡1)4e1+(24sin(1)+16cos(1)¡27)y2(x¡1)2e1+(¡2cos(1)+5=2¡2sin(1))y4e1+¡13415cos(1)+19415sin(1)¡24415¢(x¡1)5e1+¡1643¡28cos(1)¡1403sin(1)¢y2(x¡1)3e126第3章微积分问题的计算机求解+(14sin(1)¡16+10cos(1))y4(x¡1)e117试求下面级数的前n项及无穷项的和。①11£6+16£11+¢¢¢+1(5n¡4)(5n+1)+¢¢¢②µ12+13¶+µ122+132¶+¢¢¢+µ12n+13n¶+¢¢¢【求解】下面的语句可以直接求解级数的和。>>symsnk;symsum(1/(5*k-4)/(5*k+1),k,1,n)ans=-1/5/(5*n+1)+1/5>>symsum(1/(5*k-4)/(5*k+1),k,1,inf)ans=1/5>>symsnk;symsum(1/2^k+1/3^k,k,1,n)ans=-2*(1/2)^(n+1)-3/2*(1/3)^(n+1)+3/2>>symsum(1/2^k+1/3^k,k,1,inf)ans=3/2当然，无穷级数的和还可以通过极限的方式求出。18试求出下面的极限。①limn!1·122¡1+142¡1+162¡1+¢¢¢+1(2n)2¡1¸②limn!1nµ1n2+p+1n2+2p+1n2+3p+¢¢¢+1n2+np¶【求解】①可以用下面两种方法求解。>>symskn;symsum(1/((2*k)^2-1),k,1,inf)ans=1/2>>limit(symsum(1/((2*k)^2-1),k,1,n),n,inf)ans=1/2②可以由下面的语句直接求解。>>symsknlimit(n*symsum(1/(n^2+k*pi),k,1,n),n,inf)ans=第3章微积分问题的计算机求解27119试对下面数值描述的函数求取各阶数值微分，并用梯形法求取定积分。xi00.10.20.30.40.50.60.70.80.911.11.2yi02.20773.20583.44353.2412.81642.3111.81011.36020.981720.679070.44730.27684【求解】可以由下面的语句得出函数的各阶导数，得出的曲线如图3-2所示。>>x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2];y=[0,2.2077,3.2058,3.4435,3.241,2.8164,2.311,1.8101,...1.3602,0.9817,0.6791,0.4473,0.2768];[dy1,dx1]=diff_ctr(y,x(2)-x(1),1);[dy2,dx2]=diff_ctr(y,x(2)-x(1),2);[dy3,dx3]=diff_ctr(y,x(2)-x(1),3);[dy4,dx4]=diff_ctr(y,x(2)-x(1),4);plot(dx1+x(1),dy1,'-',dx2+x(1),dy2,'--',dx3+x(1),dy3,':',dx4+x(1),dy4,'-.')0.30.40.50.60.70.80.911.11.2−200−150−100−50050100n=1n=2图3-2各阶导数的数值解曲线20试求出以下的曲线积分。①Zl(x2+y2)ds，l为曲线x=a(cost+tsint);y=a(sint¡tcost);(06t62p)。②Zl(yx3+ey)dx+(xy3+xey¡2y)dy，其中l为a2x2+b2y2=c2正向上半椭圆。③Zlydx¡xdy+(x2+y2)dz，l为曲线x=et;y=e¡t;z=at，06t61，a>0。④Zl(exsiny¡my)dx+(excosy¡m)dy，其中l为由(a;0)点到(0;0)再经x2+y2=ax上正向半圆周构成的曲线。【求解】套用书中给出的第一类和第二类曲线积分公式，则可以直接得出曲线积分的结果。28第3章微积分问题的计算机求解>>symsat;x=a*(cos(t)+t*sin(t));y=a*(sin(t)-t*cos(t));f=x^2+y^2;I=int(f*sqrt(diff(x,t)^2+diff(y,t)^2),t,0,2*pi)I=2*a^3*pi^2+4*a^3*pi^4>>symsxyabct;x=c*cos(t)/a;y=c*sin(t)/b;P=y*x^3+exp(y);Q=x*y^3+x*exp(y)-2*y;ds=[diff(x,t);diff(y,t)];I=int([PQ]*ds,t,0,pi)I=2/15*c*(2*c^4-15*b^4)/a/b^4>>symst;symsapositive;x=exp(t);y=exp(-t);z=a*t;F=[y,-x,(x^2+y^2)];ds=[diff(x,t);diff(y,t);diff(z,t)];I=int(F*ds,t,0,1)I=2+1/2*a*exp(1)^2-1/2*a*exp(-1)^2>>symstm;symsapositive;x1=t;y1=0;F1=[exp(x1)*sin(y1)-m*y1,exp(x1)*cos(y1)-m];x2=a/2+a/2*cos(t);y2=a/2*sin(t);F2=[exp(x2)*sin(y2)-m*y2,exp(x2)*cos(y2)-m];I1=int(F1*[diff(x1,t);diff(y1,t)],t,0,a)I2=int(F2*[diff(x2,t);diff(y2,t)],t,0,pi);I=I1+I2I2=1/8*a^2*m*pi21试求出下面的曲面积分。①ZSµ2x+4y3+z¶ds，S为平面x2+y3+z4=1与x=¡1;y=¡2;z=¡3围成的表面。②ZSx2y2zdxdy，其中S为半球面z=pR2¡x2¡y2的下侧。【求解】第4章线性代数问题的计算机求解1Jordan矩阵是矩阵分析中一类很实用的矩阵，其一般形式为J=26664¡®10¢¢¢00¡®1¢¢¢0...............000¢¢¢¡®37775，例如J1=266664¡510000¡510000¡510000¡510000¡5377775试利用diag()函数给出构造J1的语句。【求解】利用diag()能够构造对角矩阵和次对角矩阵的性质，可以由下面语句建立起所需的矩阵。>>J1=diag([-5-5-5-5-5])+diag([1111],1)J1=-510000-510000-510000-510000-52幂零矩阵是一类特殊的矩阵，其基本形式为Hn=2666664010¢¢¢0001¢¢¢0...............000¢¢¢1000¢¢¢03777775亦即，矩阵的次主对角线元素为1，其余均为0，试验证对指定阶次的幂零矩阵，有Hin=0对所有的i>n成立。【求解】可以用循环的方式构造出各阶幂零矩阵，并对其求出i+1次方，判定得出矩阵的范数，若发现范数大于0的阶次，则显示其阶次，若为零矩阵则不显示任何内容。通过运行下面的语句，可见不显示任何内容，故i<100的幂零矩阵满足上述性质。>>fori=1:100A=diag(ones(1,i),1);ifnorm(A^(1+i))>0,disp(i);endend30第4章线性代数问题的计算机求解3试从矩阵的显示格式区分符号矩阵和数值矩阵，明确它们的含义和应用场合。若A矩阵为数值矩阵，B为符号矩阵，C=A*B运算得出的C矩阵是符号矩阵还是数值矩阵？【求解】得出的结果当然是符号矩阵，见下例。>>A=ones(5);B=sym(A);A*Bans=[5,5,5,5,5][5,5,5,5,5][5,5,5,5,5][5,5,5,5,5][5,5,5,5,5]4请将下面给出的矩阵A和B输入到MATLAB环境中，并将它们转换成符号矩阵。A=2666666664576516523100146420644396366210760077724407748672173777777775;B=266666666435501233254625121134635152124101201¡3¡4¡7378121¡107¡68153777777775【求解】矩阵的输入与转换是很直接的。>>A=[5,7,6,5,1,6,5;2,3,1,0,0,1,4;6,4,2,0,6,4,4;3,9,6,3,6,6,2;10,7,6,0,0,7,7;7,2,4,4,0,7,7;4,8,6,7,2,1,7];A=sym(A)A=[5,7,6,5,1,6,5][2,3,1,0,0,1,4][6,4,2,0,6,4,4][3,9,6,3,6,6,2][10,7,6,0,0,7,7][7,2,4,4,0,7,7][4,8,6,7,2,1,7]>>B=[3,5,5,0,1,2,3;3,2,5,4,6,2,5;1,2,1,1,3,4,6;3,5,1,5,2,1,2;4,1,0,1,2,0,1;-3,-4,-7,3,7,8,12;1,-10,7,-6,8,1,5];B=sym(B)B=[3,5,5,0,1,2,3][3,2,5,4,6,2,5][1,2,1,1,3,4,6][3,5,1,5,2,1,2][4,1,0,1,2,0,1]第4章线性代数问题的计算机求解31[-3,-4,-7,3,7,8,12][1,-10,7,-6,8,1,5]5试求出Vandermonde矩阵A=266664a4a3b4b3b2b1c4c3d4d3d2d1e4e3e2e1377775的行列式，并以最简的形式显示结果。【求解】利用书中编写的面向符号矩阵的vander()函数，可以构造出Vandermonde矩阵并可以求出该矩阵的行列式。>>symsabcde;A=vander([abcde])A=[a^4,a^3,a^2,a,1][b^4,b^3,b^2,b,1][c^4,c^3,c^2,c,1][d^4,d^3,d^2,d,1][e^4,e^3,e^2,e,1]>>det(A),simple(ans)ans=(c-d)*(b-d)*(b-c)*(a-d)*(a-c)*(a-b)*(-d+e)*(e-c)*(e-b)*(e-a)6利用MATLAB语言提供的现成函数对习题4中给出的两个矩阵进行分析，判定它们是否为奇异矩阵，得出矩阵的秩、行列式、迹和逆矩阵，检验得出的逆矩阵是否正确。【求解】以A矩阵为例，可以对其进行如下分析。>>A=[5,7,6,5,1,6,5;2,3,1,0,0,1,4;6,4,2,0,6,4,4;3,9,6,3,6,6,2;10,7,6,0,0,7,7;7,2,4,4,0,7,7;4,8,6,7,2,1,7];A=sym(A);rank(A)ans=7>>det(A)ans=-35432>>trace(A)ans=27>>B=inv(A);latex(B)32第4章线性代数问题的计算机求解A¡1=266666666666666422974429¡1157442927718858¡170944291734429¡12784429¡17844296465885815378858415117716¡13324429¡4694429¡17584429¡14678858¡2404717716¡465117716¡206413543260918858307988584827885864391771655158858¡9778858418517716¡14244429¡7944429¡8424429¡4918858¡616317716¡887177164633543216138858¡127885812718858156717716378117716256517716219354321898858¡92188584298858¡335317716¡922517716495917716¡66633543221378858¡98858260188581785177163777777777777775>>A*Bans=[1,0,0,0,0,0,0][0,1,0,0,0,0,0][0,0,1,0,0,0,0][0,0,0,1,0,0,0][0,0,0,0,1,0,0][0,0,0,0,0,1,0][0,0,0,0,0,0,1]7试求出习题4中给出的A和B矩阵的特征多项式、特征值与特征向量，并验证Hamilton-Caylay定理，解释并验证如何运算能消除误差。【求解】仍以A矩阵为例。>>A=[5,7,6,5,1,6,5;2,3,1,0,0,1,4;6,4,2,0,6,4,4;3,9,6,3,6,6,2;10,7,6,0,0,7,7;7,2,4,4,0,7,7;4,8,6,7,2,1,7];A=sym(A);eig(A)ans=5.009396680079366526215873006955228.679593193974410579078264020229.27480714110743938760483528351799e-1+1.1755376247101009492093136044131*i-1.6336795424500642956747726147329+6.9740721596526560301948635104611*i-3.4765922173751363914655588544224-1.6336795424500642956747726147329-6.9740721596526560301948635104611*i.27480714110743938760483528351799e-1-1.1755376247101009492093136044131*i>>p=poly(A)p=x^7-27*x^6-18*x^5-1000*x^4+3018*x^3+24129*x^2+2731*x+35432>>p=sym2poly(p)p=第4章线性代数问题的计算机求解331-27-18-1000301824129273135432>>polyvalm(p,A)ans=[0,0,0,0,0,0,0][0,0,0,0,0,0,0][0,0,0,0,0,0,0][0,0,0,0,0,0,0][0,0,0,0,0,0,0][0,0,0,0,0,0,0][0,0,0,0,0,0,0]8试对习题4中给出的A和B矩阵进行奇异值分解、LU分解及正交分解矩阵。【求解】仍以A矩阵为例，这里只考虑数值解。>>A=[5,7,6,5,1,6,5;2,3,1,0,0,1,4;6,4,2,0,6,4,4;3,9,6,3,6,6,2;10,7,6,0,0,7,7;7,2,4,4,0,7,7;4,8,6,7,2,1,7];[s,v,d]=svd(A)s=-0.4239-0.10610.2030-0.4091-0.07460.05760.7691-0.14810.08940.08360.4027-0.20520.87030.0378-0.29740.0356-0.57760.48400.4