高等应用数学问题的MATLAB求解习题参考解答.pdf

第 2 章MATLAB 语言程序设计基础 1 启动 MATLAB 环境，并给出语句 tic,? A=rand(500);? B=inv(A);? norm(A*B-eye(500), toc，试运行该语句，观察得出的结果，并利用 help 命令对你不熟悉的语句进行帮助信息查 询，逐条给出上述程序段与结果的解释。 ? tic,? A=rand(500);? B=inv(A);? norm(A*B-eye(500),? toc ans? = 1.2333e-012 Elapsed? time? is? 1.301000? seconds. 2 试用符号元素工具箱支持的方式表达多项式 f(x)?= x5+?3 x4+?4 x3+? 2 x 2 +?3 x +?6 ，并令 x = s 1 s +?1 ，将 f(x) 替换成 s 的函数。 【求解】 可以先定义出 f 函数，则由 subs() 函数将 x 替换成 s 的函数 ? syms? s? x f=x5+3*x4+4*x3+2*x2+3*x+6; F=subs(f,x,(s-1)/(s+1) 3 用 MATLAB 语句输入矩阵 A 和 B 矩阵 前面给出的是 4 4 矩阵，如果给出 A(5,6)?=?5命令将得出什么结果？ 【求解】 用课程介绍的方法可以直接输入这两个矩阵 ? A=1? 2? 3? 4;? 4? 3? 2? 1;? 2? 3? 4? 1;? 3? 2? 4? 1 ? 若给出 A(5,6)=5 命令，虽然这时的行和列数均大于 B 矩阵当前的维数，但仍然可以执行该 语句，得出 ? A(5,6)=5 复数矩阵也可以用直观的语句输入 ? B=1+4i? 2+3i? 3+2i? 4+1i;?4+1i? 3+2i? 2+3i? 1+4i; 2+3i? 3+2i? 4+1i? 1+4i;? 3+2i? 2+3i? 4+1i?1+4i; 4假设已知矩阵 A，试给出相应的 MATLAB 命令，将其全部偶数行提取出来，赋给 B 矩阵， 用 A =magic(8) 命令生成 A 矩阵，用上述的命令检验一下结果是不是正确。 【求解】 魔方矩阵可以采用 magic() 生成，子矩阵也可以提取出来 ? A=magic(8),? B=A(2:2:end,:) 5 用 MATLAB 语言实现下面的分段函数 y = f(x)?= h,x?D h/Dx,| x |6 D h,x? y=h*(xD)? +? h/D*x.*(abs(x)? for? i=1:length(x) if? x(i)D,? y(i)=h; elseif? abs(x(i) ? sum(sym(2).1:63) ? 由于结果有 19 位数值，所以用 double 型不能精确表示结果，该数据类型最多表示 16 位有效 数字。其实用符号运算方式可以任意保留有效数字，例如可以求 200 项的和或 1000 项的和可 以由下面语句立即得出。 ? sum(sym(2).1:200) ? ? sum(sym(2).1:1000) 9 已知 Fibonacci 数列由式 ak= ak1+ak2, 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-函数。 function y=fib(n) if round(n)=n else, y=1; end else error(n must be positive integer.) end 例如，n = 10 可以求出相应的项为 fib(10) 现在需要比较一下递归实现的速度和循环实现的速度 tic, fib(20), toc tic, a=1 1; for i=3:30, a(i)=a(i-1)+a(i-2); end, a(30), toc 应该指出，递归的调用方式速度较慢，比循环语句慢很多，所以不是特别需要，解这样问题 没有必要用递归调用的方式。 7 编写一个矩阵相加函数 mat add()，使其具体的调用格式为 A=mat add(A1,A2,A3,)， 要求该函数能接受任意多个矩阵进行解法运算。 【求解】 可以编写下面的函数，用 varargin 变量来表示可变输入变量 function A=mat_add(varargin) A=0; for i=1:length(varargin), A=A+varargini; end 如果想得到合适的错误显示，则可以试用 try, catch 结构。 function A=mat_add(varargin) try A=0; for i=1:length(varargin), A=A+varargini; end catch, error(lasterr); end 8 自己编写一个 MATLAB 函数，使它能自动生成一个 m m 的 Hankel 矩阵，并使其调用格 式为 v=h1,h2,hm,hm+1, ,h2m1; H=myhankel(v)。 【求解】 解决这样的问题可以有多种方法： 最直接的方法，Hi,j= hi+j1，利用双重循环 function H=myhankel(v) m=(length(v)+1)/2;% 严格说来还应该判定给定输入向量长度奇偶性 for i=1:m, for j=1:m H(i,j)=v(i+j-1); end,end 考虑某一行 (或列)，ai= hi,hi+1, ,hi+m1，就可以用单重循环生成 Hankel 矩阵了 function H=myhankel(v) m=(length(v)+1)/2;% 严格说来还应该判定给定输入向量长度奇偶性 for i=1:m, H(i,:)=v(i:i+m-1); end 利用现有的 hankel() 函数，则 function H=myhankel(v) m=(length(v)+1)/2;% 严格说来还应该判定给定输入向量长度奇偶性 H=hankel(v(1:m),v(m:end); 10 由矩阵理论可知，如果一个矩阵 M 可以写成 M = A + BCBT, 并且其中 A, B, C 为相应 阶数的矩阵，则 M 矩阵的逆矩阵可以由下面的算法求出 M1= A + BCBT 1 = A1 A1B C1+ BTA1B 1 BTA1 试根据上面的算法用 MATLAB 语句编写一个函数对矩阵 M 进行求逆，并通过一个小例子来 检验该程序，并和直接求逆方法进行精度上的比较。 【求解】 编写这个函数 function Minv=part_inv(A,B,C) Minv=inv(A)-inv(A)*B*inv(inv(C)+B*inv(A)*B)*B*inv(A); 假设矩阵为 M = 51503616 50776032 36608748 16324868 且已知该矩阵可以分解成 A = 1000 0200 0030 0004 , B = 1234 2340 3400 4000 , C = 4000 0300 0020 0001 对这个例子。可以 M=51 50 36 16; 50 77 60 32; 36 60 87 48; 16 32 48 68; iM=inv(M); % 数值逆，直接解法 A=diag(1 2 3 4); B=hankel(1 2 3 4); C=diag(4 3 2 1); iM1=part_inv(A,B,C)% 分块矩阵的求解方法 乍看结果，似乎二者完全一致，实际上数值算法是有区别的。我们这里用解析方法得出矩阵 的逆，然后用下面的语句比较两个结果的精度 M1=sym(M);iM0=inv(M1) norm(double(iM0)-iM) % 直接求解的误差范数 norm(double(iM0)-iM1) % 间接求解的误差范数 可见，用间接方法得出的逆矩阵误差更大，因为在这里新编写的函数中 inv() 函数使用 了多次，由此产生很大的传递误差。由此可以得出结论：如果某问题存在直接解，则尽量别 使用间接方法，以加大传递误差。 11 下面给出了一个迭代模型 ( xk+1= 1 + yk 1.4x2 k yk+1= 0.3xk 写出求解该模型的 M-函数，如果取迭代初值为 x0= 0,y0= 0，那么请进行 30000 次迭代求出 一组 x 和 y 向量，然后在所有的 xk和 yk坐标处点亮一个点 (注意不要连线)，最后绘制出所 需的图形。提示 这样绘制出的图形又称为 Henon 引力线图，它将迭代出来的随机点吸引到一 起，最后得出貌似连贯的引力线图。 【求解】 用循环形式解决此问题，可以得出如图 2-1 所示的 Henon 引力线图。 x=0; y=0; for i=1:29999 x(i+1)=1+y(i)-1.4*x(i)2; y(i+1)=0.3*x(i); end plot(x,y,.) 上述的算法由于动态定义 x 和 y，所以每循环一步需要重新定维，这样做是很消耗时间 的，所以为加快速度，可以考虑预先定义这两个变量，如给出 x=zeros(1,30000)。 12 用 MATLAB 语言的基本语句显然可以立即绘制一个正三角形，试结合循环结构，编写一个 小程序，在同一个坐标系下绘制出该正三角形绕其中心旋转后得出的一系列三角形，还可以 调整旋转步距观察效果。 【求解】 假设正三角形逆时针旋转 度，则可以得出如图 2-2a 所示的示意图，三角形的三个 顶点为 (cos,sin), (cos( +120),sin( +120), (cos( +240),sin( +240)，可以绘制出 其曲线，如图 2-2b 所示，试减小步距，如选择 = 2,1,0.1，观察效果。 t=0,120,240,0*pi/180;% 变换成弧度 xxx=; yyy=; for i=0:5:360 tt=i*pi/180; xxx=xxx; cos(tt+t); yyy=yyy; sin(tt+t); end plot(xxx,yyy,r), axis(square) 13 选择合适的步距绘制出下面的图形 sin ?1 t ，其中 t (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) 14 对合适的 范围选取分别绘制出下列极坐标图形 = 1.00132， = 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 所示。 ezplot(x2+y2-3*x*y2); hold on ezplot(x3-x2=y2-y) 可用局部放大的方法求出更精确的值，如图 2-5b 所示。从图上可以精确读出两个交 点，(0.4012,0.8916)，(1.5894,0.8185)。试将这两个点分别代入原始方程进行验证。 16 请分别绘制出 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) 17 在图形绘制语句中，若函数值为不定式 NaN，则相应的部分不绘制出来，试利用该规律绘制 z = sinxy 的表面图，并剪切下 x2+ y26 0.52的部分。 【求解】 给出下面命令可以得出矩形区域的函数值，再找出 x2+ y26 0.52区域的坐标，将其 函数值设置成 NaN，最终得出如图 2-7 所示的曲面。 第 3 章微积分问题的计算机求解 1 试求出如下极限。 lim x(3 x + 9x) 1 x ，lim x (x + 2)x+2(x + 3)x+3 (x + 5)2x+5 【求解】 极限问题由下面的语句可以直接求出。 syms x; f=(3x+9x)(1/x); limit(f,x,inf) syms x; f=(x+2)(x+2)*(x+3)(x+3)/(x+5)(2*x+5); limit(f,x,inf) 2 试求下面的双重极限。 lim x1 y2 x2y + xy3 (x + y)3 ，lim x0 y0 xy xy + 1 1，lim x0 y0 1 cos x2 + y2 x2 + y2ex 2+y2 。 【求解】 双重极限问题可以由下面语句直接求解。 syms x y fa=(x2*y+x*y3)/(x+y)3; limit(limit(fa,x,-1),y,2) fb=x*y/(sqrt(x*y+1)-1); limit(limit(fb,x,0),y,0) fc=(1-cos(x2+y2)*exp(x2+y2)/(x2+y2); limit(limit(fc,x,0),y,0) 3 求出下面函数的导数。 y(x) = q xsinx1 ex， y = 1 cosax x(1 cos ax) atan y x = ln(x2+ y2)， y(x) = 1 na ln xn+ a xn , n 0 【求解】 由求导函数 diff() 可以直接得出如下结果，其中 为隐函数，故需要用隐函数求 导公式得出导数。 syms x; f=sqrt(x*sin(x)*sqrt(1-exp(x); simple(diff(f) syms a x y=(1-sqrt(cos(a*x)/(x*(1-cos(sqrt(a*x) simple(diff(y) f=atan(y/x)-log(x2+y2); f1=simple(-diff(f,x)/diff(f,y) f1 = (y+2*x)/(x-2*y) syms n positive; syms a; f=-log(xn+a)/xn)/(n*a); diff(f,x) 4 试求出 y(t) = s (x 1)(x 2) (x 3)(x 4) 函数的 4 阶导数。 【求解】 高阶导数可以由下面语句直接得出 syms a x f=sqrt(x-1)*(x-2)/(x-3)/(x-4); simple(diff(f,x,4) 5 在高等数学中，求解分子和分母均同时为 0 或 时，分式极限时可使用 LH opital 法则，即 对分子分母分别求导数，再由比值得出，试用该法则 lim x0 ln(1 + x)ln(1 x) ln(1 x2) x4 ， 并和直接求出的极限结果相比较。 【求解】 从给出的分母看，若想使之在 x = 0 处的值不为 0，则应该对其求 4 阶导数，同样， 还应该对分子求 4 阶导数，将 x = 0 代入结果，这样就可以使用 LH opital 法则求出极限了。 syms x; n=log(1+x)*log(1-x)-log(1-x2); d=x4; n4=diff(n,x,4); d4=diff(d,x,4); n4=subs(n4,x,0); L=n4/d4 现在直接求极限可以验证上述结果是正确的。 limit(n/d,x,0) 6 已知参数方程 x = lncost y = cost tsint ，试求出 dy dx 和 d2y dx2 fl fl fl fl t=?/3 。 【求解】 参数方程的导数可以由下面语句直接求出。 syms t; x=log(cos(t); y=cos(t)-t*sin(t); diff(y,t)/diff(x,t) f=diff(y,t,2)/diff(x,t,2); subs(f,t,sym(pi)/3) 7 假设 u = cos1 r x y ，试验证 2u xy = 2u yx。 【求解】 证明二者相等亦可以由二者之差为零来证明，故由下面的语句直接证明。 syms x y; u=acos(x/y); diff(diff(u,x),y)-diff(diff(u,y),x) 8 设 ( xu + yv = 0 yu + xv = 1 ，试求解 2u xy 。 【求解】 用下面的语句可以直接得出如下结果。 syms x y u v u,v=solve(x*u+y*v=0,y*u+x*v=1,u,v); diff(diff(u,x),y) 9 假设 f(x,y) = Z xy 0 et 2dt，试求x y 2f x2 2 2f xy + 2f y2 。 【求解】 由下面的命令可以得出所需结果。 syms x y t f=int(exp(-t2),t,0,x*y); x/y*diff(f,x,2)-2*diff(diff(f,x),y)+diff(f,y,2) simple(ans) 10 假设已知函数矩阵 f(x,y,z) = 3x + eyz x3+ y2sinz ，试求出其 Jacobi 矩阵。 【求解】 Jacobi 矩阵可以由下面的语句直接得出。 syms x y z F=3*x+exp(y)*z; x3+y2*sin(z); jacobian(F,x,y,z) 11 试求解下面的不定积分问题。 I(x) = Z 3x2+ a x2(x2+ a)2 dx， I(x) = Z px(x + 1) x +1 + xdx I(x) = Z xeaxcosbxdx， I(t) = Z eaxsinbxsincxdx 【求解】 该不定积分可以由下面的命令直接求出 syms x a f=(3*x2+a)/(x2+(x2+a)2); int(f,x) 可以用下面的语句求出问题的解 syms x; f=sqrt(x*(x+1)/(sqrt(x)+sqrt(x+1); int(f,x); latex(ans) 可以求出下面的结果 syms a b x f=x*exp(a*x)*cos(b*x); int(f,x); latex(ans) 用下面的语句求解，得 syms x a b c; f=exp(a*x)*sin(b*x)*sin(c*x); latex(int(f,x) 12 试求出下面的定积分或无穷积分。 I = Z 0 cosx xdx， I = Z 1 0 1 + x2 1 + x4 dx 【求解】 可以直接求解 syms x; int(cos(x)/sqrt(x),x,0,inf) 可以得出 syms x; int(1+x2)/(1+x4),x,0,1) 13 假设 f(x) = e5xsin(3x + ?/3)，试求出积分函数 R(t) = Z t 0 f(x)f(t + x)dx。 【求解】 定义了 x 的函数，则可以由 subs() 函数定义出 t + x 的函数，这样由下面的语句可 以直接得出 R 函数。 syms x t; f=exp(-5*x)*sin(3*x+sym(pi)/3); R=int(f*subs(f,x,t+x),x,0,t); simple(R) 14 对 a 的不同取值试求出 I = Z 0 cosax 1 + x2 dx。 【求解】 由下面的循环结构可以得出不同 a 值下的无穷积分值，并可以绘制出它们之间关系 的曲线，如图 3-1 所示。 syms x a; f=cos(a*x)/(1+x2); aa=0:0.1:pi; y=; for n=aa b=int(subs(f,a,n),x,0,inf); y=y, double(b); end plot(aa,y) 15 试对下面函数进行 Fourier 幂级数展开。 f(x) = (? |x|)sinx, ? 6 x syms x; f=(sym(pi)-abs(x)*sin(x); A,B,F=fseries(f,x,10,-pi,pi); F 可以由下面语句求解，并得出数学公式为 syms x; f=exp(abs(x); A,B,F=fseries(f,x,10,-pi,pi); F 似乎求解起来更困难，巧妙利用符号运算工具箱中的 heaviside() 函数，则可以将原函数 表示成 f(x) = 2 heaviside x ? 2 2 ?x |x ?/2| x ?/2 这样就可以用下面的语句求出函数的 Fourier 级数。 syms x; 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); F 16 试求出下面函数的 Taylor 幂级数展开。 Z x 0 sint t dt ln ?1 + x 1 x ln x + p 1 + x2 (1 + 4.2x2)0.2 e5xsin(3x + ?/3) 分别关于 x = 0、x = a 的幂级数展开 对 f(x,y) = 1 cos x2 + y2 x2 + y2ex 2+y2 关于 x = 1,y = 0 进行二维 Taylor 幂级数展开。 【求解】 由下面的语句可以分别求出各个函数的幂级数展开，由 latex(ans) 函数可以得出下 面的数学表示形式。 syms t x; f=int(sin(t)/t,t,0,x); taylor(f,x,15) syms x; f=log(1+x)/(1-x), taylor(f,x,15) syms x; f=log(x+sqrt(1+x2); taylor(f,x,15) syms x; f=(1+4.2*x2)0.2; taylor(f,x,13) syms x a; f=exp(-5*x)*sin(3*x+sym(pi)/3); taylor(f,x,4,a) 该函数需要使用 Maple 的展开函数。 syms x y; f=(1-cos(x2+y2)/(x2+y2)*exp(x2+y2); F=maple(mtaylor,f,x=1,y,4) 17 试求下面级数的前 n 项及无穷项的和。 1 1 6 + 1 6 11 + + 1 (5n 4)(5n + 1) + ?1 2 + 1 3 + ? 1 22 + 1 32 + + ? 1 2n + 1 3n + 【求解】 下面的语句可以直接求解级数的和。 syms n k; symsum(1/(5*k-4)/(5*k+1),k,1,n) symsum(1/(5*k-4)/(5*k+1),k,1,inf) syms n k; symsum(1/2k+1/3k,k,1,n) symsum(1/2k+1/3k,k,1,inf) 18 试求出下面的极限。 lim n 1 22 1 + 1 42 1 + 1 62 1 + + 1 (2n)2 1 lim nn ? 1 n2+ ? + 1 n2+ 2? + 1 n2+ 3? + + 1 n2+ n? 【求解】 可以用下面两种方法求解。 syms k n; symsum(1/(2*k)2-1),k,1,inf) limit( symsum(1/(2*k)2-1),k,1,n),n,inf) 可以由下面的语句直接求解。 syms k n limit(n*symsum(1/(n2+k*pi),k,1,n),n,inf) 19 试对下面数值描述的函数求取各阶数值微分，并用梯形法求取定积分。 xi00.10.20.30.40.50.60.70.80.911.11.2 yi02.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,-.) 第 4 章线性代数问题的计算机求解 1 Jordan 矩阵是矩阵分析中一类很实用的矩阵，其一般形式为 J = 100 010 . . . . . . . . . . . 000 ，例如 J1= 51000 05100 00510 00051 00005 试利用 diag() 函数给出构造 J1的语句。 【求解】 利用 diag() 能够构造对角矩阵和次对角矩阵的性质，可以由下面语句建立起所需的 矩阵。 J1=diag(-5 -5 -5 -5 -5)+diag(1 1 1 1,1) 2 幂零矩阵是一类特殊的矩阵，其基本形式为 Hn= 0100 0010 . . . . . . . . . . . 0001 0000 亦即，矩阵的次主对角线元素为 1，其余均为 0，试验证对指定阶次的幂零矩阵，有 Hi n= 0 对所有的 i n 成立。 【求解】 可以用循环的方式构造出各阶幂零矩阵，并对其求出 i + 1 次方，判定得出矩阵的范 数，若发现范数大于 0 的阶次，则显示其阶次，若为零矩阵则不显示任何内容。通过运行下 面的语句，可见不显示任何内容，故 i for i=1:100 A=diag(ones(1,i),1); if norm(A(1+i)0, disp(i); end end 3 试从矩阵的显示格式区分符号矩阵和数值矩阵，明确它们的含义和应用场合。若 A 矩阵为数 值矩阵，B 为符号矩阵，C=A*B 运算得出的 C 矩阵是符号矩阵还是数值矩阵？ 【求解】 得出的结果当然是符号矩阵，见下例。 A=ones(5); B=sym(A); A*B 4 请将下面给出的矩阵 A 和 B 输入到 MATLAB 环境中，并将它们转换成符号矩阵。 A = 5765165 2310014 6420644 3963662 10760077 7244077 4867217 ,B = 3550123 3254625 1211346 3515212 4101201 34737812 11076815 【求解】 矩阵的输入与转换是很直接的。 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) 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) 5 试求出 Vandermonde 矩阵 A = a4a3a2a1 b4b3b2b1 c4c3c2c1 d4d3d2d1 e4e3e2e1 的行列式，并以最简的形式显示结果。 【求解】 利用书中编写的面向符号矩阵的 vander() 函数，可以构造出 Vandermonde 矩阵并 可以求出该矩阵的行列式。 syms a b c d e; A=vander(a b c d e) det(A), simple(ans) 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) det(A) trace(A) B=inv(A); latex(B) A*B 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) p=poly(A) p=sym2poly(p) polyvalm(p,A) 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) L U=lu(A) orth(A) 9 试求出下面矩阵的特征值、特征向量、奇异值。 A = 27577 74933 39838 59636 26854 , B = 703795980137661 54795727112284 445523252894469 69588087619965 621173737299988 【求解】 A 矩阵的分析如下。 A=2,7,5,7,7; 7,4,9,3,3; 3,9,8,3,8; 5,9,6,3,6; 2,6,8,5,4 B=703,795,980,137,661; 547,957,271,12,284; 445,523,252,894,469; B=16,17,9,12,12; 17,12,12,2,18; 9,12,18,7,13; 12,2,7,18,12; 12,18,13,12,10; chol(B) 11 试对矩阵 A = 20.50.50.5 01.50.50.5 20.54.50.5 2122 进行 Jordan 变换，并得出变换矩阵。 【求解】 将 A 输入，并转换成符号矩阵，这样就能得出 Jordan 变换矩阵与 Jordan 矩阵。 A=-2,0.5,-0.5,0.5; 0,-1.5,0.5,-0.5; 2,0.5,-4.5,0.5; 2,1,-2,-2; V J=jordan(sym(A) 12 试求下面齐次方程的基础解系。 6x1+ x2+ 4x3 7x4 3x5= 0 2x1 7x2 8x3+ 6x4= 0 4x1+ 5x2+ x3 6x4+ 8x5= 0 34x1+ 36x2+ 9x3 21x4+ 49x5= 0 26x1 12x2 27x3+ 27x4+ 17x5= 0 【求解】 可以将方程写成矩阵形式，得出的两列向量为方程的基础解系。 A=6,1,4,-7,-3; -2,-7,-8,6,0; -4,5,1,-6,8; -34,36,9,-21,49; -26,-12,-27,27,17; A=sym(A); rank(A) null(A) 13 试求下面线性代数方程的解析解与数值解，并检验解的正确性。 29321 1011050 82463 56684 X = 140 384 033 953 【求解】 求出 A, A,B 两个矩阵的秩，可见二者相同，所以方程不是矛盾方程，应该有无穷 多解。 A=2,-9,3,-2,-1; 10,-1,10,5,0; 8,-2,-4,-6,3; -5,-6,-6,-8,-4; B=-1,-4,0; -3,-8,-4; 0,3,3; 9,-5,3; rank(A), rank(A B) 用下面的语句可以求出方程的解析解，并可以验证该解没有误差。 x0=null(sym(A); x_analytical=sym(A)B; syms a; x=a*x0 x0 x0+x_analytical A*x-B 用数值解方法也可以求出方程的解，但会存在误差，且与任意常数 a 的值有关。 x0=null(A); x_numerical=AB; syms a; x=a*x0 x0 x0+x_numerical; vpa(x,10) A*x-B 14 试判定下面的线性代数方程是否有解。 162313 511108 97612 414151 X = 1 3 4 7 【求解】 由秩判定矩阵可以得出如下结果。 A=16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1; B=1; 3; 4; 7; rank(A), rank(A B) 15 试求出线性代数方程的解析解，并验证解的正确性。 29412586 12287337 30357510 31166991 11214687 5181911118 262710151318 X = 19 512 412 109 05 1018 202 【求解】 先判定方程解的存在性，得出结论，方程有无穷多解。 A=2,9,4,12,5,8,6; 12,2,8,7,3,3,7; 3,0,3,5,7,5,10; 3,11,6,6,9,9,1; 11,2,1,4,6,8,7; 5,-18,1,-9,11,-1,18; 26,-27,-1,0,-15,-13,18; X=1,9; 5,12; 4,12; 10,9; 0,5; 10,18; -20,2; rank(A), rank(A X) 这样就能解出方程的解析解，并验证方程解是这正确的。 x0=null(sym(A); syms a1 a2; x1=sym(A)X; x=a1*x0(:,1) x0(:,1)+a2*x0(:,2) x0(:,2)+x1 A*x-X 16 试用数值方法和解析方法求取下面的 Sylvester 方程，并验证得出的结果。 36405 14224 63673 13100110 04034 X + X 321 292 219 = 211 412 561 644 663 【求解】 用下面的语句可以直接求出其数值解，并可以判定该解的精度。 A=3,-6,-4,0,5; 1,4,2,-2,4; -6,3,-6,7,3; -13,10,0,-11,0; 0,4,0,3,4; B=3,-2,1; -2,-9,2; -2,-1,9; C=-2,1,-1; 4,1,2; 5,-6,1; 6,-4,-4; -6,6,-3; X=lyap(A,B,C) norm(A*X+X*B+C) 利用书中给出的解析解方法，可以求出该方程的解析解。 X=lyap(sym(A),B,C) A*X+X*B+C 18 假设已知某 Jordan 块矩阵 A 及其组成部分为 A= A1 A2 A3 , A1= 310 031 003 , A2= 5 1 05 , A3= 1100 0110 0011 0001 试用解析解运算的方式得出 eAt，sin 2At + ? 3 ，eA 2tA2 + sin(A3t)At + esinAt。 【求解】 用下面的语句可以求出矩阵问题的解析解。 A1=diag(-3 -3 -3)+diag(1 1,1); A2=-5 1; 0 -5; A3=diag(-1 -1 -1 -1)+diag(1 1 1,1); A=sym(diagm(A1,A2,A3); syms t; expm(A*t) syms x t; funm(A,sin(2*x*t+sym(pi)/3),x) funm(A,exp(x2*t)*x2+sin(x3*t)*x*t+exp(sin(x*t),x) 例如可以求出上述结果的第一行第二列元素为 54te9t 6e9t 81cos(27t)t2 sin(27t)t + cos(3t)tesin(3t) 19 假设已知矩阵 A 如下给出，试求出 eAt,sinAt,eAtsin A2eAtt。 A = 4.500.51.5 0.540.50.5 1.512.51.5 0113 【求解】 用下面的语句可以直接得出矩阵指数解析解。 A=-4.5,0,0.5,-1.5; -0.5,-4,0.5,-0.5; 1.5,1,-2.5,1.5; 0,-1,-1,-3; A=sym(A); syms t; expm(A*t) syms x t; funm(A,sin(x*t),x) syms x t; funm(A,exp(x*t)*sin(x2*exp(x*t)*t),x) 第 5 章积分变换与复变函数问题的计算机求解 1 对下列的函数 f(t) 进行 Laplace 变换 fa(t) = sint t , fb(t) = t5sint, fc(t) = t8cost fd(t) = t6et fe(t) = 5eat+ t4eat+ 8e2t, ff(t) = etsin(t + ), fg(t) = e12t+ 6e9t 【求解】 给定函数的 Laplace 变换可以用下面的语句直接得出 syms a t; f=sin(a*t)/t; laplace(f) syms t a; f=t5*sin(a*t); laplace(f) syms t a; f=t8*cos(a*t); laplace(f) syms t a; f=t6*exp(a*t); laplace(f) syms a t; f=5*exp(-a*t)+t4*exp(-a*t)+8*exp(-2*t); laplace(f) syms a b c t; f=exp(b*t)*sin(a*t+c); laplace(f) syms t; f=exp(-12*t)+6*exp(9*t); laplace(f) 2 对下面的 F(s) 式进行 Laplace 反变换 Fa(s) = 1 s2(s2 a2)(s + b) , Fb(s) = s a s b, F c(s) = ln s a s b Fd(s) = 1 s(