Matlab 参考教程第五章 数值计算文档格式.docx
《Matlab 参考教程第五章 数值计算文档格式.docx》由会员分享，可在线阅读，更多相关《Matlab 参考教程第五章 数值计算文档格式.docx（75页珍藏版）》请在冰点文库上搜索。
eri=
3.0708e-004
rei=
6.6280e-007
（3）“左除”法解恰定方程的误差、残差、运算次数和所用时间
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
实际误差er_actual'
）,disp（er_actual）,disp（'
'
）
近似的最大可能误差er_approx'
）,disp（er_approx）,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的比较。
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;
%使中间迭代结果不显示。
[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
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.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
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
Columns1through4
-0.0328+0.0243i000
0-0.0328-0.0243i00
000.0130+0.0379i0
0000.0130-0.0379i
0000
Column5
0
0.0396
1.0e+010*
2.1016
2.0251
1.9796
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_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
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
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（:
21）,Z（:
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
Dx_2=
11
FX=
111
FY=
FX_2=
222
FY_2=
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,'
-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'
%画坐标横轴
y（t）'
图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
y4=
exitflag=
1
（5）求在tt（3）附近的精确零点
从理论分析可知，
是函数的一个零点。
但即便是以十分靠近该零点的
为搜索的初始值，也找不到
，而却找到了另一个零点。
原因是曲线没有穿越横轴。
请看下面指令的运行结果。
[t3,y3,exitflag]=fzero（y,tt（3）,[],a,b）
[0.58266,-0.59706].
t3=
-0.5198
y3=
（6）观察fzero所采用的options缺省设置，并更改控
升级会员 