高等应用数学问题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-1 MATLAB 演示程序界面
1
例如，用户选择MATLAB ! Graphics ! Volume Vlsulization 演示，则将得出如图1-2 所示的
演示说明，单击其中的Run this demo 栏目，则将得出如图1-3 所示的演示界面。用户可以在
该界面下按按钮，逐步演示相关内容，而实现这样演示的语句将在该程序界面的下部窗口中 给出。
2 作者用MATLAB 语言编写了给出例子的源程序，读者可以自己用type 语句阅读一下源程
序，对照数学问题初步理解语句的含义，编写的源程序说明由下表列出。
第1 章计算机数学语言概述3 图1-2 MATLAB 演示程序界面举例 序号文件名程序说明
例1.1 c1ex1.m 利用MATLAB 的符号运算工具箱求解微分问题 例1.2 c1ex2.m 分别利用MATLAB 的符号运算工具箱和数值运算功能求解多项式方程，其中用数值方法得出 的结果有误差
例1.3 c1ex3.m 分别利用MATLAB 的符号运算工具箱和数值运算功能计算Hilbert 矩阵的行列式，其中用数值 方法得出的结果有很大误差
例1.4 c1ex4.m 令x1 = y; x2 = y_，则可以将原来的二阶微分方程转换成
2
一阶微分方程组，然后就可以求解微分
方程的数值解了，原方程是非线性微分方程，故不存在解析解。ode45() 函数可以求解常微分方
程组，而dde23() 可以求解延迟微分方程，或更直观地采用Simulink 绘制求解框图。
例1.5 c1ex5.m 线性规划问题调用最优化工具箱中的linprog() 函数可以立即得出结果，若想求解整数规划问 题，则需要首先安装整数规划程序ipslv mex()。 4 第1 章计算机数学语言概述 图1-3 MATLAB 体视化演示程序界面 第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)), toc
3
ans = 1.2333e-012
Elapsed time is 1.301000 seconds.
2 试用符号元素工具箱支持的方式表达多项式f(x) = x5 + 3x4 + 4x3 + 2x2 + 3x + 6，并令
x = s ? 1 s + 1
，将f(x) 替换成s 的函数。
【求解】可以先定义出f 函数，则由subs() 函数将x 替换成s 的函数 >> syms s x
f=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)+6 3 用MATLAB 语句输入矩阵A 和B 矩阵 ① A = 2 664 1 2 3 4 4 3 2 1 2 3 4 1
4
3 2 4 1 3 775; ②
B = 2 664
1 + 4j 2 + 3j 3 + 2j 4 + 1j 4 + 1j 3 + 2j 2 + 3j 1 + 4j 2 + 3j 3 + 2j 4 + 1j 1 + 4j 3 + 2j 2 + 3j 4 + 1j 1 + 4j 3 775
前面给出的是4 ￡ 4 矩阵，如果给出A(5; 6) = 5 命令将得出什么结果？
【求解】用课程介绍的方法可以直接输入这两个矩阵 >> A=[1 2 3 4; 4 3 2 1; 2 3 4 1; 3 2 4 1] A = 1 2 3 4
6 第2 章MATLAB 语言程序设计基础 4 3 2 1
5
2 3 4 1 3 2 4 1
若给出A(5,6)=5 命令，虽然这时的行和列数均大于B 矩阵当前的维数，但仍然可以执行该 语句，得出 >> A(5,6)=5 A = 1 2 3 4 0 0 4 3 2 1 0 0 2 3 4 1 0 0 3 2 4 1 0 0 0 0 0 0 0 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]; B =
1.0000 + 4.0000i 2.0000 + 3.0000i 3.0000 + 2.0000i 4.0000 + 1.0000i 4.0000 + 1.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 1.0000 + 4.0000i 2.0000 + 3.0000i 3.0000 + 2.0000i 4.0000 + 1.0000i 1.0000 + 4.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 4.0000 + 1.0000i 1.0000 + 4.0000i 4 假设已知矩阵A，试给出相应的MATLAB 命令，将其全部偶数行提取出来，赋给B 矩阵，
6
用A =magic(8) 命令生成A 矩阵，用上述的命令检验一下结果是不是正确。
【求解】魔方矩阵可以采用magic() 生成，子矩阵也可以提取出来 >> A=magic(8), B=A(2:2:end,:) A =
64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1 B =
9 55 54 12 13 51 50 16
第2 章MATLAB 语言程序设计基础7 40 26 27 37 36 30 31 33 41 23 22 44 45 19 18 48 8 58 59 5 4 62 63 1
5 用MATLAB 语言实现下面的分段函数y = f(x) = 8< :
7
h; x > D h=Dx; j x j6 D ?h; x < ?D 。
【求解】两种方法，其一，巧用比较表达式解决 >> y=h*(x>D) + h/D*x.*(abs(x)<=D) -h*(x<-D); 另外一种方法，用循环语句和条件转移语句 >> for i=1:length(x) if x(i)>D, y(i)=h;
elseif abs(x(i))<=D, y(i)= h/D*x(i); else, y(i)=-h; end end
其中，前者语句结构简单，但适用范围更广，允许使用矩阵型x，后者只能使用向量型的
x，但不能处理矩阵问题。
6 用数值方法可以求出S = X63
i=0
2i = 1 + 2 + 4 + 8 + ￠ ￠ ￠ + 262 + 263，试不采用循环的形式求出 和式的数值解。由于数值方法采用double 形式进行计算的，难以保证有效位数字，所以结果
不一定精确。试采用符号运算的方法求该和式的精确值。
8
【求解】用符号运算的方式可以采用下面语句 >> 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 =
214301721437253464189685009812000362112280962341106721488750077674070
210224987224498639675763139171625518934583510629365037429057138462808
8 第2 章MATLAB 语言程序设计基础
719691551493971496078691355496484619708421492101247422837559083643060
929499671638825347975351183310878921541258291423929553730843353
9
208596
63305248773674411336138750
7 编写一个矩阵相加函数mat add()，使其具体的调用格式为A=mat add(A1,A2,A3,￠ ￠ ￠ )，
要求该函数能接受任意多个矩阵进行解法运算。
【求解】可以编写下面的函数，用varargin 变量来表示可变输入变量 function A=mat_add(varargin) A=0;
for i=1:length(varargin), A=A+varargin{i}; end
如果想得到合适的错误显示，则可以试用try, catch 结构。 function A=mat_add(varargin) try A=0;
for i=1:length(varargin), A=A+varargin{i}; end catch, error(lasterr); end
8 自己编写一个MATLAB 函数，使它能自动生成一个m ￡ m 的Hankel 矩阵，并使其调用格
式为v=[h1; h2; hm; hm+1; ￠ ￠ ￠ ; h2m?1]; H=myhankel(v)。 【求解】解决这样的问题可以有多种方法： ①最直接的方法，Hi;j = hi+j?1，利用双重循环 function H=myhankel(v)
m=(length(v)+1)/2; % 严格说来还应该判定给定输入向量长度奇偶性
10
for i=1:m, for j=1:m H(i,j)=v(i+j-1); end,end
②考虑某一行(或列)，ai = [hi; hi+1; ￠ ￠ ￠ ; hi+m?1]，就可以用单重循环生成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));
第2 章MATLAB 语言程序设计基础9
9 已知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 项，这可以使用递归调用的功能，
11
而递归调用的出口为
1。综上，可以编写出M-函数。 function y=fib(n) if round(n)==n & n>=1 if n>=3
y=fib(n-1)+fib(n-2); else, y=1; end else
error('n must be positive integer.') end
例如，n = 10 可以求出相应的项为 >> fib(10) ans = 55
现在需要比较一下递归实现的速度和循环实现的速度 >> tic, fib(20), toc ans = 832040 elapsed_time = 62.0490
>> tic, a=[1 1]; for i=3:30, a(i)=a(i-1)+a(i-2); end, a(30), toc ans =
12
832040 elapsed_time = 0.0100
应该指出，递归的调用方式速度较慢，比循环语句慢很多，所以不是特别需要，解这样问题 没有必要用递归调用的方式。 10 第2 章MATLAB 语言程序设计基础
10 由矩阵理论可知，如果一个矩阵M 可以写成M = A + BCBT , 并且其中A, B, C 为相应
阶数的矩阵，则M 矩阵的逆矩阵可以由下面的算法求出
M?1 =
3
A + BCBT ′?1 = A?1 ? A?1B 3
C?1 + BTA?1B ′?1
BTA?1
试根据上面的算法用MATLAB 语句编写一个函数对矩阵M 进行求逆，并通过一个小例子来
检验该程序，并和直接求逆方法进行精度上的比较。
13
【求解】编写这个函数 function Minv=part_inv(A,B,C)
Minv=inv(A)-inv(A)*B*inv(inv(C)+B'*inv(A)*B)*B'*inv(A); 假设矩阵为
M =
2 664 51 50 36 16 50 77 60 32 36 60 87 48 16 32 48 68 3 775
且已知该矩阵可以分解成
A =
2 664 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 3
14
775
; B =
2 664 1 2 3 4 2 3 4 0 3 4 0 0 4 0 0 0 3 775
; C =
2 664 4 0 0 0 0 3 0 0 0 0 2 0 0 0 0 1 3 775
对这个例子。可以
>> M=[51 50 36 16; 50 77 60 32; 36 60 87 48; 16 32 48 68]; iM=inv(M); % 数值逆，直接解法
15
iM =
0.0553 -0.0389 0.0017 0.0041 -0.0389 0.0555 -0.0210 -0.0021 0.0017 -0.0210 0.0328 -0.0137 0.0041 -0.0021 -0.0137 0.0244
>> A=diag([1 2 3 4]); B=hankel([1 2 3 4]); C=diag([4 3 2 1]); iM1=part_inv(A,B,C) % 分块矩阵的求解方法 iM1 =
0.0553 -0.0389 0.0017 0.0041 -0.0389 0.0555 -0.0210 -0.0021 0.0017 -0.0210 0.0328 -0.0137 0.0041 -0.0021 -0.0137 0.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]
16
>> norm(double(iM0)-iM) % 直接求解的误差范数 ans = 2.7990e-017
>> norm(double(iM0)-iM1) % 间接求解的误差范数 ans = 3.6583e-016
可见，用间接方法得出的逆矩阵误差更大，因为在这里新编写的函数中inv() 函数使用
了多次，由此产生很大的传递误差。由此可以得出结论：如果某问题存在直接解，则尽量别
使用间接方法，以加大传递误差。 11 下面给出了一个迭代模型(
xk+1 = 1 + yk ? 1:4x2 k yk+1 = 0:3xk 写出求解该模型的M-函数，如果取迭代初值为x0 = 0,y0 = 0，那么请进行30000 次迭代求出
一组x 和y 向量，然后在所有的xk 和yk 坐标处点亮一个点(注意不要连线)，最后绘制出所
需的图形。提示这样绘制出的图形又称为Henon 引力线图，它将迭代出来的随机点吸引到一
起，最后得出貌似连贯的引力线图。
17
【求解】用循环形式解决此问题，可以得出如图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 语言的基本语句显然可以立即绘制一个正三角形，试结合循环结构，编写一个
小程序，在同一个坐标系下绘制出该正三角形绕其中心旋转后得出的一系列三角形，还可以 调整旋转步距观察效果。
12 第2 章MATLAB 语言程序设计基础 ?1.5 ?1 ?0.5 0 0.5 1 1.5 ?0.4 ?0.3 ?0.2
18
?0.1 0 0.1 0.2 0.3 0.4
图2-1 Henon 引力线图
【求解】假设正三角形逆时针旋转μ 度，则可以得出如图2-2a 所示的示意图，三角形的三个
顶点为(cos μ; sin μ), (cos(μ +120〒); sin(μ +120〒)), (cos(μ +240〒); sin(μ +240〒))，可以绘制出
其曲线，如图2-2b 所示，试减小步距，如选择￠μ = 2; 1; 0:1，观察效果。
.............................................................
μ x y - 6
(a) 示意图 ?1 ?0.5 0 0.5 1 ?1 ?0.5
19
0 0.5 1
(b) 曲线绘制效果 图2-2 曲线绘制
>> 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
第2 章MATLAB 语言程序设计基础13 plot(xxx',yyy','r'), axis('square')
13 选择合适的步距绘制出下面的图形sin μ 1
t ?
，其中t 2 (?1; 1)。
【求解】用普通的绘图形式，选择等间距，得出如图2-3a 所示的曲线，其中x = 0 左右显得 粗糙。
20
a2 + (b + c)2
12 试求出下面的定积分或无穷积分。 ① I = Z 1 0 cos x p x dx， ② I = Z 1 0 1 + x2 1 + x4 dx 【求解】① 可以直接求解 >> syms x; int(cos(x)/sqrt(x),x,0,inf) ans =
1/2*2^(1/2)*pi^(1/2) ② 可以得出
>> syms x; int((1+x^2)/(1+x^4),x,0,1) ans = 1/4*2^(1/2)*pi
13 假设f(x) = e?5x sin(3x + p=3)，试求出积分函数R(t) =
46
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) 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 该结果可以写成 1 1360 15 (et)10 p 3 cos (3t) ? 68 cos (3t) ? 15 (et)10 sin (3t) ? 25
p 3 sin (9t) + 25 (et)10 p 3 sin (3t)
+15 sin (9t) ? 25 cos (9t) ? 15
47
p 3 cos (9t) + 93 (et)10 cos (3t) (et)15
14 对a 的不同取值试求出I = Z 1 0 cos ax 1 + x2 dx。
【求解】由下面的循环结构可以得出不同a 值下的无穷积分值，并可以绘制出它们之间关系 的曲线，如图3-1 所示。
第3 章微积分问题的计算机求解23
>> syms x a; f=cos(a*x)/(1+x^2); 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)
0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6
48
0.8 1 1.2 1.4 1.6
图3-1 不同a 值下的积分值曲线
15 试对下面函数进行Fourier 幂级数展开。
① f(x) = (p ? jxj) sin x; ?p 6 x < p ② f(x) = ejxj; ?p 6 x < p ③ f(x) = (
2x=l; 0 < x < l=2 2(l ? x)=l; l=2 < x < l ，且l = p。
【求解】① 可以立即由下面的语句求出。 >> syms x; f=(sym(pi)-abs(x))*sin(x); [A,B,F]=fseries(f,x,10,-pi,pi); F F =
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 下可以显示为 1 2p sin x +
49
16 9 sin 2x p + 32 225 sin 4x p + 48 1225 sin 6x p + 64 3969 sin 8x p + 80 9801
50