数值计算方法精讲.docx
- 文档编号:12061004
- 上传时间:2023-04-16
- 格式:DOCX
- 页数:87
- 大小:602.06KB
数值计算方法精讲.docx
《数值计算方法精讲.docx》由会员分享,可在线阅读,更多相关《数值计算方法精讲.docx(87页珍藏版)》请在冰豆网上搜索。
数值计算方法精讲
数值计算方法
5.1引言
本章将花较大的篇幅讨论若干常见数值计算问题:
线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。
但与一般数值计算教科书不同,本章的讨论重点是:
如何利用现有的世界顶级数值计算资源MATLAB。
至于数学描述,本章将遵循“最低限度自封闭”的原则处理,以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。
对于那些熟悉其他高级语言(如FORTRAN,Pascal,C++)的读者来说,通过本章,MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗,感受到摆脱烦琐编程后的眉眼舒展。
对于那些经过大学基本数学教程的读者来说,通过本章,MATLAB精良完善的计算指令,自然易读的程序将使他们感悟“教程”数学的基础地位和局限性,看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。
对于那些熟悉MATLAB基本指令的读者来说,通过本章,围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系,获得综合运用不同指令解决具体问题的思路和借鉴。
由于MATLAB的基本运算单元是数组,所以本章内容将从矩阵分析、线性代数的数值计算开始。
然后再介绍函数零点、极值的求取,数值微积分,数理统计和分析,拟合和插值,Fourier分析,和一般常微分方程初值问题。
本章的最后讨论稀疏矩阵的处理,因为这只有在大型问题中,才须特别处理。
从总体上讲,本章各节之间没有依从关系,即读者没有必要从头到尾系统阅读本章内容。
读者完全可以根据需要阅读有关节次。
除特别说明外,每节中的例题指令是独立完整的,因此读者可以很容易地在自己机器上实践。
5.1LU分解和恰定方程组的解
5.1.1LU分解、行列式和逆
(1)LU分解
(2)行列式和逆
5.1.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.1.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.2矩阵特征值和矩阵函数
5.2.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.7098e-010
5.2.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.2.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.2.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.3奇异值分解
5.3.1奇异值分解和矩阵结构
5.3.1.1奇异值分解定义
5.3.1.2矩阵结构的奇异值分解描述
(1)秩
(2)零空间(Nullspace)和值空间(Rangespace)
(3)范数(Norm)
(4)矩阵条件数
(5)子空间夹角
(6)广义逆(Moore-Penrosepseudoinverse)
5.3.2线性二乘问题的解
5.3.2.1矩阵除运算的广义化
5.3.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.4函数的数值导数和切平面
5.4.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.4.2偏导数和梯度
5.4.2.1理论定义
5.4.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.5函数的零点
5.5.1多项式的根
5.5.2一元函数的零点
5.5.2.1利用MATLAB作图指令获取初步近似解
5.5.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
y4=
0
exitflag=
1
(5)求在tt(3)附近的精确零点
从理论分析可知,
是函数的一个零点。
但即便是以十分靠近该零点的
为搜索的初始值,也找不到
,而却找到了另一个零点。
原因是曲线没有穿越横轴。
请看下面指令的运行结果。
[t3,y3,exitflag]=fzero(y,tt(3),[],a,b)
Zerofoundintheinterval:
[0.58266,-0.59706].
t3=
-0.5198
y3=
0
exitflag=
1
(6)观察fzero所采用的options缺省设置,并更改控制计算精度的相对误差设置。
op=optimset('fzero')%提取fz
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法