04数值计算.docx

04数值计算.docx

《04数值计算.docx》由会员分享，可在线阅读，更多相关《04数值计算.docx（75页珍藏版）》请在冰点文库上搜索。

04数值计算

第四章数值计算

4.1引言

本章将花较大的篇幅讨论若干常见数值计算问题：

线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。

但与一般数值计算教科书不同，本章的讨论重点是：

如何利用现有的世界顶级数值计算资源MATLAB。

至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。

对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。

对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。

对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。

由于MATLAB的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。

然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值、边值问题。

本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。

从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。

读者完全可以根据需要阅读有关节次。

除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。

MATLAB从5.3版升级到6.x版后，本章内容的变化如下：

●MATLAB从6.0版起，其矩阵和特征值计算指令不再以LINPACK和EISPACK库为基础，而建筑在计算速度更快、运行更可靠的LAPACK和ARPACK程序库的新基础上。

因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。

这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。

对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。

本书新版全部算例结果是在6.x版上给出的。

●在5.3版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。

这种调用方式在6.x版中已被宣布为“过渡期内允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。

因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。

●MATLAB从6.0版起，提供了一组专门求微分方程“边值问题”数值解的指令。

适应这种变化，本章新增第4.14.5节，用2个算例阐述求解细节。

●5.3版中的积分指令quad8已经废止；6.x版启用新积分指令quadl；6.5版新增三重积分指令triplequad。

本章新版对此作了相应的改变。

4.2LU分解和恰定方程组的解

4.2.1LU分解、行列式和逆

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=[68101214];

fork=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）

formatshorte

disp（'实际误差er_actual'）,disp（er_actual）,disp（''）

disp（'近似的最大可能误差er_approx'）,disp（er_approx）,disp（''）

disp（'最大可能误差er_max'）,disp（er_max）,disp（''）

Hilbert矩阵阶数

68101214

实际误差er_actual

1.5410e-0101.7310e-0071.9489e-0049.1251e-0022.1257e+000

近似的最大可能误差er_approx

3.3198e-0093.3879e-0063.5583e-0033.9846e+0009.0475e+001

最大可能误差er_max

7.9498e-0073.8709e-0021.2703e+0034.7791e+0074.0622e+010

4.3矩阵特征值和矩阵函数

4.3.1特征值和特征向量的求取

【例4.3.1-1】简单实阵的特征值问题。

A=[1,-3;2,2/3];[V,D]=eig（A）

V=

0.77460.7746

0.0430-0.6310i0.0430+0.6310i

D=

0.8333+2.4438i0

00.8333-2.4438i

【例4.3.1-2】本例演示：

如矩阵中有元素与截断误差相当时的特性值问题。

A=[3-2-0.92*eps

-24-1-eps

-eps/4eps/2-10

-0.5-0.50.11];

[V1,D1]=eig（A）;ER1=A*V1-V1*D1

[V2,D2]=eig（A,'nobalance'）;ER2=A*V2-V2*D2

ER1=

0.00000.00000.00000.0000

0-0.0000-0.0000-0.0000

0.0000-0.0000-0.00000.0000

0.00000.00000.0000-0.5216

ER2=

1.0e-014*

-0.26650.0111-0.0559-0.1055

0.44410.12210.03430.0833

0.00220.00020.00070

0.0194-0.02220.02220.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

formatshorte,c_equ=cond（B）,c_eig=condeig（B）

B=

1000

0100

0011

0001

c_equ=

2.6180e+000

Warning:

Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=1.110223e-016.

>InD:

\MATLAB6P1\toolbox\matlab\matfun\condeig.matline30

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=

-911-2163-252

70-69141-4211684

-575575-11493451-13801

3891-38917782-2334593365

1024-10242048-614424572

DJ=

01000

00100

00010

00001

00000

D=

Columns1through4

-0.0408000

0-0.0119+0.0386i00

00-0.0119-0.0386i0

0000.0323+0.0230i

0000

Column5

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.77460

0.0430-0.6310

DR=

0.83332.4438

-2.44380.8333

4.3.4矩阵的谱分解和矩阵函数

【例4.3.4-1】数组乘方与矩阵乘方的比较。

clear,A=[123;456;789];

A_Ap=A.^0.3

A_Mp=A^0.3

A_Ap=

1.00001.23111.3904

1.51571.62071.7118

1.79281.86611.9332

A_Mp=

0.6962+0.6032i0.4358+0.1636i0.1755-0.2759i

0.6325+0.0666i0.7309+0.0181i0.8292-0.0305i

0.5688-0.4700i1.0259-0.1275i1.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.30000.09000.0270

0.00810.00240.0007

0.00020.00010.0000

pA_M=

2.93420.4175-1.0993

-0.02780.7495-0.4731

-1.9898-0.91841.1531

【例4.3.4-3】sin的数组运算和矩阵运算比较。

（A取自例4.3.4-1）

A_sinA=sin（A）

A_sinM=funm（A,'sin'）

A_sinA=

0.84150.90930.1411

-0.7568-0.9589-0.2794

0.65700.98940.4121

A_sinM=

-0.6928-0.23060.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.77.56.30.83-0.082]';

x=inv（A'*A）*A'*y,xx=pinv（A）*y,xxx=A\y

Warning:

Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.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:

Rankdeficient,rank=3tol=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=

123

456

789

Dx=

333

333

Dx_2=

11

11

11

FX=

111

111

111

FY=

333

333

333

FX_2=

222

222

222

FY_2=

666

666

666

【例4.5.2.2-2】研究偶极子（Dipole）的电势（Electricpotential）和电场强度（Electricfielddensity）。

设在

处有电荷

，在

处有电荷

。

那么在电荷所在平面上任何一点的电势和场强分别为

，

。

其中

。

。

又设电荷

，

，

。

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}偶极子的场'）,holdon

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'）,holdoff

图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'）;holdon,plot（t,zeros（size（t））,'k'）;

xlabel（'t'）;ylabel（'y（t）'）,holdoff

图4.6-1

（3）

由于Notebook中无法实现zoom、ginput指令涉及的图形和鼠标交互操作，因此下面指令必须在MATLAB指令窗中运行，并得到如图4.6-2所示的局部放大图及鼠标操作线。

zoomon

[tt,yy]=ginput（5）;zoomoff

图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）

holdon,contour（X,Y,F2,v）,holdoff

图4.6-3

（2）

[x0,y0]=ginput

（2）;

disp（[x0,y0]）

-0.7926-0.7843

0.79260.7843

（3）

fun='[sin（x

（1）-x

（2））,cos（x

（1）+x

（2））]';%<12>

[xy,f,exit]=fsolve（fun,[x0

（2）,y0

（2）]）%<13>

Optimizationterminatedsuccessfully:

First-orderoptimalitylessthanOPTIONS.TolFun,andnonegative/zerocurvaturedetected

xy=

0.78540.7854

f=

1.0e-006*

-0.09840.2011

exit=

1

〖说明〗

[fun.m]

functionff=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.00001.0000

sfval=

8.1