第四章数值计算.doc

第四章 数值计算

4.1 引言
本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。
对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。
对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。
对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。
由于MATLAB的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。
从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。
MATLAB从5.3版升级到6.x版后，本章内容的变化如下：
l MATLAB从6.0版起，其矩阵和特征值计算指令不再以LINPACK和EISPACK库为基础，而建筑在计算速度更快、运行更可靠的LAPACK和ARPACK程序库的新基础上。因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在6.x版上给出的。
l 在5.3版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用方式在6.x版中已被宣布为“过渡期内允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。
l MATLAB从6.0版起，提供了一组专门求微分方程“边值问题”数值解的指令。适应这种变化，本章新增第4.14.5节，用2个算例阐述求解细节。
l 5.3版中的积分指令quad8已经废止；6.x版启用新积分指令quadl ；6.5版新增三重积分指令triplequad。本章新版对此作了相应的改变。
4.2 LU分解和恰定方程组的解
4.2.1 LU分解、行列式和逆
4.2.2 恰定方程组的解
【例4.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比
（1）
randn('state',0);
A=gallery('randsvd',100,2e13,2);
x=ones(100,1);
b=A*x;
cond(A)
ans =
1.9990e+013

（2）
tic
xi=inv(A)*b;
ti=toc
eri=norm(x-xi)
rei=norm(A*xi-b)/norm(b)
ti =
0.7700
eri =
0.0469
rei =
0.0047

（3）
tic;xd=A\b;
td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)
td =
0
erd =
0.0078
red =
2.6829e-015

4.2.3 范数、条件数和方程解的精度
【例4.2.3-1】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确逆矩阵的生成函数。本例将对方程近似解和准确解进行比较。
N=[6 8 10 12 14];
for k=1:length(N)
n=N(k);
H=hilb(n);
Hi=invhilb(n);
b=ones(n,1);
x_approx=H\b;
x_exact=Hi*b;
ndb=norm(H*x_approx-b);nb=norm(b);
ndx=norm(x_approx - x_exact);nx=norm(x_approx);
er_actual(k)=ndx/nx;
K=cond(H);
er_approx(k)=K*eps;
er_max(k)=K*ndb/nb;
end

disp('Hilbert矩阵阶数'),disp(N)
format short e
disp('实际误差 er_actual'),disp(er_actual),disp('')
disp('近似的最大可能误差 er_approx'),disp(er_approx),disp('')
disp('最大可能误差 er_max'),disp(er_max),disp('')
Hilbert矩阵阶数
6 8 10 12 14
实际误差 er_actual
1.5410e-010 1.7310e-007 1.9489e-004 9.1251e-002 2.1257e+000
近似的最大可能误差 er_approx
3.3198e-009 3.3879e-006 3.5583e-003 3.9846e+000 9.0475e+001
最大可能误差 er_max
7.9498e-007 3.8709e-002 1.2703e+003 4.7791e+007 4.0622e+010

4.3 矩阵特征值和矩阵函数
4.3.1 特征值和特征向量的求取
【例4.3.1-1】简单实阵的特征值问题。
A=[1,-3;2,2/3];[V,D]=eig(A)
V =
0.7746 0.7746
0.0430 - 0.6310i 0.0430 + 0.6310i
D =
0.8333 + 2.4438i 0
0 0.8333 - 2.4438i

【例4.3.1-2】本例演示：如矩阵中有元素与截断误差相当时的特性值问题。
A=[3 -2 -0.9 2*eps
-2 4 -1 -eps
-eps/4 eps/2 -1 0
-0.5 -0.5 0.1 1 ];
[V1,D1]=eig(A);ER1=A*V1-V1*D1
[V2,D2]=eig(A,'nobalance');ER2=A*V2-V2*D2
ER1 =
0.0000 0.0000 0.0000 0.0000
0 -0.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000 0.0000
0.0000 0.0000 0.0000 -0.5216
ER2 =
1.0e-014 *
-0.2665 0.0111 -0.0559 -0.1055
0.4441 0.1221 0.0343 0.0833
0.0022 0.0002 0.0007 0
0.0194 -0.0222 0.0222 0.0333

【例4.3.1-3】指令eig与eigs的比较。
rand('state',1),A=rand(100,100)-0.5;
t0=clock;[V,D]=eig(A);T_full=etime(clock,t0)
options.tol=1e-8;
options.disp=0;
t0=clock;[v,d]=eigs(A,1,'lr',options);
T_part=etime(clock,t0)
[Dmr,k]=max(real(diag(D)));
d,D(1,1)
T_full =
0.2200
T_part =
3.1300
d =
3.0140 + 0.2555i
ans =
3.0140 + 0.2555i

vk1=V(:,k);
vk1=vk1/norm(vk1);v=v/norm(v);
V_err=acos(norm(vk1'*v))*180/pi
D_err=abs(D(k,k)-d)/abs(d)
V_err =
1.2074e-006
D_err =
4.2324e-010

4.3.2 特征值问题的条件数
【例4.3.2-1】矩阵的代数方程条件数和特征值条件数。
B=eye(4,4);B(3,4)=1;B
format short e,c_equ=cond(B),c_eig=condeig(B)
B =
1 0 0 0
0 1 0 0
0 0 1 1
0 0 0 1
c_equ =
2.6180e+000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.110223e-016.
> In D:\MATLAB6P1\toolbox\matlab\matfun\condeig.m at line 30
c_eig =
1.0000e+000
1.0000e+000
4.5036e+015
4.5036e+015

【例4.3.2-2】对亏损矩阵进行Jordan分解。
A=gallery(5)
[VJ,DJ]=jordan(A);
[V,D,c_eig]=condeig(A);c_equ=cond(A);
DJ,D,c_eig,c_equ
A =
-9 11 -21 63 -252
70 -69 141 -421 1684
-575 575 -1149 3451 -13801
3891 -3891 7782 -23345 93365
1024 -1024 2048 -6144 24572
DJ =
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
0 0 0 0 0
D =
Columns 1 through 4
-0.0408 0 0 0
0 -0.0119 + 0.0386i 0 0
0 0 -0.0119 - 0.0386i 0
0 0 0 0.0323 + 0.0230i
0 0 0 0
Column 5
0
0
0
0
0.0323 - 0.0230i
c_eig =
1.0e+010 *
2.1293
2.0796
2.0796
2.0020
2.0020
c_equ =
2.0253e+018

4.3.3 复数特征值对角阵与实数块特征值对角阵的转化
【例4.3.3-1】把例4.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。
A=[1,-3;2,2/3];[V,D]=eig(A);
[VR,DR]=cdf2rdf(V,D)
VR =
0.7746 0
0.0430 -0.6310
DR =
0.8333 2.4438
-2.4438 0.8333

4.3.4 矩阵的谱分解和矩阵函数
【例4.3.4-1】数组乘方与矩阵乘方的比较。
clear,A=[1 2 3;4 5 6;7 8 9];
A_Ap=A.^0.3
A_Mp=A^0.3
A_Ap =
1.0000 1.2311 1.3904
1.5157 1.6207 1.7118
1.7928 1.8661 1.9332
A_Mp =
0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i
0.6325 + 0.0666i 0.7309 + 0.0181i 0.8292 - 0.0305i
0.5688 - 0.4700i 1.0259 - 0.1275i 1.4830 + 0.2150i

【例4.3.4- 2】标量的数组乘方和矩阵乘方的比较。（A取自例4.3.4-1）
pA_A=(0.3).^A
pA_M=(0.3)^A
pA_A =
0.3000 0.0900 0.0270
0.0081 0.0024 0.0007
0.0002 0.0001 0.0000
pA_M =
2.9342 0.4175 -1.0993
-0.0278 0.7495 -0.4731
-1.9898 -0.9184 1.1531

【例4.3.4-3】sin的数组运算和矩阵运算比较。（A取自例4.3.4-1）
A_sinA=sin(A)
A_sinM=funm(A,'sin')
A_sinA =
0.8415 0.9093 0.1411
-0.7568 -0.9589 -0.2794
0.6570 0.9894 0.4121
A_sinM =
-0.6928 -0.2306 0.2316
-0.1724 -0.1434 -0.1143
0.3479 -0.0561 -0.4602

4.4 奇异值分解
4.4.1 奇异值分解和矩阵结构
4.4.1.1 奇异值分解定义
4.4.1.2 矩阵结构的奇异值分解描述
4.4.2 线性二乘问题的解
4.4.2.1 矩阵除运算的广义化
4.4.2.2 线性模型的最小二乘解
【例4.4.2.2-1】对于超定方程，进行三种解法比较。其中取MATLAB库中的特殊函数生成。
（1）
A=gallery(5);A(:,1)=[];y=[1.7 7.5 6.3 0.83 -0.082]';
x=inv(A'*A)*A'*y,xx=pinv(A)*y,xxx=A\y
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 5.405078e-018.
x =
3.4811e+000
5.1595e+000
9.5340e-001
-4.6569e-002
xx =
3.4759e+000
5.1948e+000
7.1207e-001
-1.1007e-001
Warning: Rank deficient, rank = 3 tol = 1.0829e-010.
xxx =
3.4605e+000
5.2987e+000
0
-2.9742e-001

（2）
nx=norm(x),nxx=norm(xx),nxxx=norm(xxx)
nx =
6.2968e+000
nxx =
6.2918e+000
nxxx =
6.3356e+000

（3）
e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx)
e =
6.9863e-001
ee =
4.7424e-002
eee =
4.7424e-002

4.5 函数的数值导数和切平面
4.5.1 法线
【例4.5.1-1】曲面法线演示。
y=-1:0.1:1;x=2*cos(asin(y));
[X,Y,Z]=cylinder(x,20);
surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21));
view([120,18])

图 4.5.1-1
4.5.2 偏导数和梯度
4.5.2.1 理论定义
4.5.2.2 数值计算指令
【例4.5.2.2-1】用一个简单矩阵表现diff和gradient指令计算方式。
F=[1,2,3;4,5,6;7,8,9]
Dx=diff(F)
Dx_2=diff(F,1,2)
[FX,FY]=gradient(F)
[FX_2,FY_2]=gradient(F,0.5)
F =
1 2 3
4 5 6
7 8 9
Dx =
3 3 3
3 3 3
Dx_2 =
1 1
1 1
1 1
FX =
1 1 1
1 1 1
1 1 1
FY =
3 3 3
3 3 3
3 3 3
FX_2 =
2 2 2
2 2 2
2 2 2
FY_2 =
6 6 6
6 6 6
6 6 6

【例 4.5.2.2-2】研究偶极子(Dipole)的电势（Electric potential）和电场强度（Electric field density）。设在处有电荷，在处有电荷。那么在电荷所在平面上任何一点的电势和场强分别为，。其中。。又设电荷，，。
clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x;
[X,Y]=meshgrid(x,y);
rp=sqrt((X-a).^2+(Y-b).^2);rm=sqrt((X+a).^2+(Y+b).^2);
V=q*k*(1./rp-1./rm);
[Ex,Ey]=gradient(-V);
AE=sqrt(Ex.^2+Ey.^2);Ex=Ex./AE;Ey=Ey./AE;
cv=linspace(min(min(V)),max(max(V)),49);
contourf(X,Y,V,cv,'k-')
%axis('square')
title('\fontname{隶书}\fontsize{22}偶极子的场'),hold on
quiver(X,Y,Ex,Ey,0.7)
plot(a,b,'wo',a,b,'w+')
plot(-a,-b,'wo',-a,-b,'w-')
xlabel('x');ylabel('y'),hold off

图 4.5.2.2-1
4.6 函数的零点
4.6.1 多项式的根
4.6.2 一元函数的零点
4.6.2.1 利用MATLAB作图指令获取初步近似解
4.6.2.2 任意一元函数零点的精确解
【例 4.6.2.2-1】通过求的零点，综合叙述相关指令的用法。
（1）
y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b'); %<1>

（2）
a=0.1;b=0.5;t=-10:0.01:10;
y_char=vectorize(y); % <3>
Y=feval(y_char,t,a,b);
clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k');
xlabel('t');ylabel('y(t)'),hold off

图4.6-1

（3）
由于Notebook中无法实现zoom、ginput指令涉及的图形和鼠标交互操作，因此下面指令必须在MATLAB指令窗中运行，并得到如图4.6-2所示的局部放大图及鼠标操作线。
zoom on
[tt,yy]=ginput(5);zoom off

图 4.6-2

tt
tt =
-2.0032
-0.5415
-0.0072
0.5876
1.6561

（4）
[t4,y4,exitflag]=fzero(y,tt(4),[],a,b) %<11>
t4 =
0.5993
y4 =
0
exitflag =
1

（5）
[t3,y3,exitflag]=fzero(y,tt(3),[],a,b)
t3 =
-0.5198
y3 =
-5.5511e-017
exitflag =
1

（6）
op=optimset('fzero')
op =
ActiveConstrTol: []
......
Display: 'notify'
......
TolX: 2.2204e-016
TypicalX: []

op=optimset('tolx',0.01);
op.TolX
ans =
0.0100

（7）
[t4n,y4n,exitflag]=fzero(y,tt(4),op,a,b)
t4n =
0.6042
y4n =
0.0017
exitflag =
1

4.6.3 多元函数的零点
【例 4.6.3-1】求解二元函数方程组的零点。
（1）
x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y);
F1=sin(X-Y);F2=cos(X+Y);
v=[-0.2, 0, 0.2];
contour(X,Y,F1,v)
hold on,contour(X,Y,F2,v),hold off

图4.6-3

（2）
[x0,y0]=ginput(2);
disp([x0,y0])
-0.7926 -0.7843
0.7926 0.7843

（3）
fun='[sin(x(1)-x(2)),cos(x(1)+x(2))]'; %<12>
[xy,f,exit]=fsolve(fun,[x0(2),y0(2)]) %<13>
Optimization terminated successfully:
First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected
xy =
0.7854 0.7854
f =
1.0e-006 *
-0.0984 0.2011
exit =
1

〖说明〗
[FUN.M]
function ff=fun(x)
ff(1)=sin(x(1)-x(2));
ff(2)=cos(x(1)+x(2));
4.7 函数极值点
4.7.1 一元函数的极小值点
4.7.2 多元函数的极小值点
【例4.7.2-1】求的极小值点。它即是著名的Rosenbrock's "Banana" 测试函数。该测试函数有一片浅谷，许多算法难以越过此谷。（演示本例搜索过程的文件名为exm04072_1_1.m 。）
（1）
ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');
（2）
x0=[-1.2,1];[sx,sfval,sexit,soutput]=fminsearch(ff,x0)
sx =
1.0000 1.0000
sfval =
8.1777e-010
sexit =
1
soutput =
iterations: 85
funcCount: 159
algorithm: 'Nelder-Mead simplex direct search'

（3）
[ux,sfval,uexit,uoutput,grid,hess]=fminunc(ff,x0)
Warning: Gradient must be provided for trust-region method;
using line-search method instead.

> In D:\MATLAB6P1\toolbox\optim\fminunc.m at line 211

Optimization terminated successfully:
Current search direction is a descent direction, and magnitude of
directional derivative in search direction less than 2*options.TolFun
ux =
1.0000 1.0000
sfval =
1.9116e-011
uexit =
1
uoutput =
iterations: 26
funcCount: 162
stepsize: 1.2992
firstorderopt: 5.0020e-004
algorithm: 'medium-scale: Quasi-Newton line search'
grid =
1.0e-003 *
-0.5002
-0.1888
hess =
820.4028 -409.5496
-409.5496 204.7720

4.8 数值积分
4.8.1 一元函数的数值积分
4.8.1.1 闭型数值积分
【例 4.8.1.1-1】求，其精确值为 。
（1）
syms x;IS=int('exp(-x*x)','x',0,1)
vpa(IS)
IS =
1/2*erf(1)*pi^(1/2)
ans =
.74682413281242702539946743613185

（2）
fun=inline('exp(-x.*x)','x');
Isim=quad(fun,0,1),IL=quadl(fun,0,1)
Isim =
0.7468
IL =
0.7468

（3）
Ig=gauss10(fun,0,1)
Ig =
0.7463

（4）
xx=0:0.1:1.5;ff=exp(-xx.^2);
pp=spline(xx,ff);
int_pp=fnint(pp);
Ssp=ppval(int_pp,[0,1])*[-1;1]
Ssp =
0.7468

（5）

图4.8-1

4.8.1.2 开型数值积分
[GAUSS10.M]
function g = gauss10(fun,a,b)
%GAUSS10(fun,a,b)
% fun
%======================================================
x = [0.1488743390;0.4333953941;0.6974095683;...
0.8650633667;0.9739065285];
w = [0.2955242247;0.2692667193;0.2190863625;...
0.1494513492;0.0666713443];
t = .5*(b+a)+.5*(b-a)*[-flipud(x);x];
W = [flipud(w);w];
g = sum(W.*feval(fun,t))*(b-a)/2;

【例 4.8.1.2-1】当时，比较解析积分和近似积分。
（1）
syms x;F=int('cos(x)','x',-1,1),vpa(F)
F =
2*sin(1)
ans =
1.6829419696157930133050046432606
（2）
aF=cos(1/sqrt(3))+cos(-1/sqrt(3))
aF =
1.6758

【例4.8.1.2-2】求，准确结果是 。
（1）
syms x;IS=vpa(int('sqrt(log(1/x))','x',0,1))
Warning: Explicit integral could not be found.
> In D:\MATLAB6P5\toolbox\symbolic\@sym\int.m at line 58
In D:\MATLAB6P5\toolbox\symbolic\@char\int.m at line 9
IS =
.88622692545275801364908374167057
（2）用quad指令求积
ff=inline('sqrt(log(1./x))','x');Isim=quad(ff,0,1)
Warning: Divide by zero.
> In D:\MATLAB6P5\toolbox\matlab\funfun\inlineeval.m at line 13
In D:\MATLAB6P5\toolbox\matlab\funfun\@inline\feval.m at line 34
In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 62
Isim =
0.8862

（3）
Ig=gauss10(ff,0,1)
Ig =
0.8861

4.8.2 多重数值积分
4.8.2.1 积分限为常数的二重积分指令

【例 4.8.2.1-1】计算和。
（1）
syms x y
ssx01=vpa(int(int(x^y,x,0,1),y,1,2))
ssx12=vpa(int(int(x^y,y,0,1),x,1,2))
Warning: Explicit integral could not be found.
> In D:\MATLAB6P5\toolbox\symbolic\@sym\int.m at line 58
ssx01 =
.40546510810816438197801311546435
ssx12 =
1.2292741343616127837489278679215

（2）
zz=inline('x.^y','x','y');
nsx01=dblquad(zz,0,1,1,2)
nsx12=dblquad(zz,1,2,0,1)
nsx01 =
0.4055
nsx12 =
1.2293

4.8.2.2 内积分限为函数的二重积分
[double_int.m]
function SS=double_int(fun,innlow,innhi,outlow,outhi)
%double_int
%
%fun
%innlow,innhi
%outlow,outhi
y1=outlow;y2=outhi;x1=innlow;x2=innhi;f_p=fun;
SS=quad(@G_yi,y1,y2,[],[],x1,x2,f_p);

[G_YI.M]
function