4 matlab字符串数组元胞数组和构架数组.docx

4 matlab字符串数组元胞数组和构架数组.docx

《4 matlab字符串数组元胞数组和构架数组.docx》由会员分享，可在线阅读，更多相关《4 matlab字符串数组元胞数组和构架数组.docx（98页珍藏版）》请在冰点文库上搜索。

4matlab字符串数组元胞数组和构架数组

5      数值计算

5.1     引言

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

线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。

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

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

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

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

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

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

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

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

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

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

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

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

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

5.2.1   LU分解、行列式和逆

（1）LU分解

（2）行列式和逆

5.2.2   恰定方程组的解

【*例5.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比

（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程。

rand（'state',12）;%选定随机种子，目的是可重复产生随机阵A。

A=rand（100,100）+1.e8;%rand（100,100）生成（100×100）均匀分布随机矩阵。

%每个随机阵元素加

的目的是使A阵条件数升高。

x=ones（100,1）;%令解向量x为全1的100元列向量。

b=A*x;%为使Ax=b方程一致，用A和x生成b向量。

cond（A）%求A阵的条件数。

ans=

1.4426e+012

（2）“求逆”法解恰定方程的误差、残差、运算次数和所用时间

flops（0）;tic%浮点运算计数器置0；启动计时器StopwatchTimer

xi=inv（A）*b;%xi是用“求逆”法解恰定方程所得的解。

ti=toc%关闭计时器，并显示解方程所用的时间。

ci=flops%“求逆”法解方程所用的运算次数

eri=norm（x-xi）%解向量xi与真解向量x的范-2误差。

rei=norm（A*xi-b）/norm（b）%方程的范-2相对残差

ti=

0.9300

ci=

2070322

eri=

3.0708e-004

rei=

6.6280e-007

（3）“左除”法解恰定方程的误差、残差、运算次数和所用时间

flops（0）;tic;xd=A\b;%是用“左除”法解恰定方程所得的解。

td=toc,cd=flops,erd=norm（x-xd）,red=norm（A*xd-b）/norm（b）

td=

0.2200

cd=

741872

erd=

3.2243e-004

red=

2.0095e-016

5.2.3   范数、条件数和方程解的精度

【*例5.2.3-1】Hilbert矩阵是著名的病态矩阵。

MATLAB中有专门的Hilbert矩阵及其准确逆矩阵的生成函数。

本例将对方程

近似解和准确解进行比较。

所谓n阶Hilbert矩阵的形式是：

。

N=[68101214];%本例计算的矩阵阶数

fork=1:

length（N）

n=N（k）;%矩阵的阶

H=hilb（n）;%产生n阶Hilbert矩阵

Hi=invhilb（n）;%产生完全准确的n阶逆Hilbert矩阵

b=ones（n,1）;%生成n阶全1向量

x_approx=H\b;%利用左除H求近似解

x_exact=Hi*b;%利用准确逆Hilbert矩阵求准确解

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）;%计算Hilbert矩阵的条件数

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

5.0339e-0118.5981e-0082.2819e-0041.3381e-0013.9641e+000

近似的最大可能误差er_approx

3.3198e-0093.3879e-0063.5583e-0033.9259e+0003.4573e+002

最大可能误差er_max

6.0095e-0072.4531e-0021.4094e+0032.9206e+0072.4178e+010

5.3     矩阵特征值和矩阵函数

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

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

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

V=

-0.7728+0.0527i-0.7728-0.0527i

0+0.6325i0-0.6325i

D=

0.8333+2.4438i0

00.8333-2.4438i

【*例5.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-0.0000-0.00000.0000

0.0000-0.00000.00000.0000

0.00000.0000-0.00000.0000

-0.00000.00000.00001.1227

ER2=

1.0e-015*

0.4441-0.22200.1471-0.2220

00.0555-0.36290.2776

-0.0172-0.00150.00660

0-0.2220-0.11100.1388

【*例5.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）%指令eig的运作时间。

options.tol=1e-8;%为eigs设定计算精度。

options.disp=0;%使中间迭代结果不显示。

t0=clock;[v,d]=eigs（A,1,'lr',options）;%计算最大实部特征值和特征向量。

T_part=etime（clock,t0）%指令eigs的运作时间。

[Dmr,k]=max（real（diag（D）））;%在eig求得的全部特征值中找最大实部的那个。

d,D（1,1）

T_full=

2.8000

T_part=

19.5500

d=

3.0140-0.2555i

ans=

3.0140+0.2555i

vk1=V（:

k+1）;%与d相同的特征向量应是V的第k+1列。

vk1=vk1/norm（vk1）;v=v/norm（v）;%向量长度归一。

V_err=acos（norm（vk1'*v））*180/pi%求复数向量之间的夹角（度）。

D_err=abs（D（k+1,k+1）-d）/abs（d）%求两个特征值间的相对误差。

V_err=

8.5377e-007

D_err=

1.1.   7098e-010

5.3.2   特征值问题的条件数

【例5.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.500000e-018.

>InE:

\MAT53\toolbox\matlab\matfun\condeig.matline30

c_eig=

1.0000e+000

1.0000e+000

3.3333e+017

3.3333e+017

【*例5.3.2-2】对亏损矩阵进行Jordan分解。

A=gallery（5）%MATLAB设置的特殊矩阵，它具有五重特征值。

[VJ,DJ]=jordan（A）;%求出准确的特征值，使A*VJ=VJ*D成立。

[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.0328+0.0243i000

0-0.0328-0.0243i00

000.0130+0.0379i0

0000.0130-0.0379i

0000

Column5

0

0

0

0

0.0396

c_eig=

1.0e+010*

2.1016

2.1016

2.0251

2.0251

1.9796

c_equ=

5.2133e+017

5.3.3   复数特征值对角阵与实数块特征值对角阵的转化

【*例5.3.3-1】把例5.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。

[VR,DR]=cdf2rdf（V,D）

VR=

-0.77280.0527

00.6325

DR=

0.83332.4438

-2.44380.8333

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

【*例5.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

【*例5.3.4-2】标量的数组乘方和矩阵乘方的比较。

（A取自例5.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

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

（A取自例5.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

5.4     奇异值分解

5.4.1   奇异值分解和矩阵结构

5.4.1.1 奇异值分解定义

5.4.1.2 矩阵结构的奇异值分解描述

（1）秩

（2）零空间（Nullspace）和值空间（Rangespace）

（3）范数（Norm）

（4）矩阵条件数

（5）子空间夹角

（6）广义逆（Moore-Penrosepseudoinverse）

5.4.2   线性二乘问题的解

5.4.2.1 矩阵除运算的广义化

5.4.2.2 线性模型的最小二乘解

【*例5.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=7.087751e-018.

x=

3.4811e+000

5.1595e+000

9.5342e-001

-4.6562e-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=

1.1020e+000

ee=

4.7424e-002

eee=

4.7424e-002

5.5     函数的数值导数和切平面

5.5.1   法线

【*例5.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]）%控制观察角

图5.5.1-1旋转半椭球面和法线

5.5.2   偏导数和梯度

5.5.2.1 理论定义

5.5.2.2 数值计算指令

【*例5.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）%相邻列差分。

第三输入宗量2表示“列”差分。

[FX,FY]=gradient（F）%数据点步长默认为1

[FX_2,FY_2]=gradient（F,0.5）%数据点步长为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

【*例5.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）;%产生49个电位值

contourf（X,Y,V,cv,'k-'）%用黑实线画填色等位线图

%axis（'square'）%在Notebook中，此指令不用

title（'\fontname{隶书}\fontsize{22}偶极子的场'）,holdon

quiver（X,Y,Ex,Ey,0.7）%第五输入宗量0.7使场强箭头长短适中。

plot（a,b,'wo',a,b,'w+'）%用白线画正电荷位置

plot（-a,-b,'wo',-a,-b,'w-'）%用白线画负电荷位置

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

图5.5.2.2-1电偶极子的场和等位线

5.6     函数的零点

5.6.1   多项式的根

5.6.2   一元函数的零点

5.6.2.1 利用MATLAB作图指令获取初步近似解

5.6.2.2 任意一元函数零点的精确解

【*例5.6.2.2-1】通过求

的零点，综合叙述相关指令的用法。

（1）构造一个内联函数对象

被解函数

以

为自变量，

和

为参数。

假如在fzero中直接采用字符串表示被解函数，容易出错。

因此先构造内联函数如下：

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）;%为避免循环，把y改写成适合数组运算形式。

<4>

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

图5.6.2.2-1函数零点分布观察图

（3）利用zoom和ginput指令获得零点的初始近似值（在MATLAB指令窗中进行）

zoomon%在MATLAB指令窗中运行，获局部放大图

[tt,yy]=ginput（5）;zoomoff%在MATLAB指令窗中运行，用鼠标获5个零点猜测值。

图5.6.2.2-2局部放大和利用鼠标取值图

tt%显示所得零点初始猜测值（该指令可在Notebook中运行）。

tt=

-2.0032

-0.5415

-0.0072

0.5876

1.6561

（4）求靠近tt（4）的精确零点

[t4,y4,exitflag]=fzero（y,tt（4）,[],a,b）%<11>

Zerofoundintheinterval:

[0.57094,0.60418].

t4=

0。

5993