现代数值计算（第3版）课件全套 第1-9章 科学计算与MATLAB-常微分方程初边值问题数值解

第一章AdvancedNumericalComputing科学计算与MATLAB现代数值计算同济大学数学科学学院目录/Contents第一章科学计算与MATLAB第一节科学计算的意义第二节误差基础知识第三节MATLAB软件1.1科学计算的意义科学计算的出现利用现代计算机辅助，解决实际问题计算的挑战：基因测序，全球天气模拟科学计算问题的主要步骤数学建模数值算法评价科学计算软件MATLAB，Mathematica,Maple,Python,1.2误差基础知识实际问题数学问题(模型误差)计算问题(截断误差、观测误差)结果(舍入误差)1.2.1误差的来源实际问题数学模型数值算法结果编程处理模型误差观测误差舍入误差截断误差比较检验1.2误差基础知识设有真值,及近似值，称为该近似值的绝对误差1.2.2误差度量称为绝对误差限若称为相对误差称为相对误差限由于真值难以求出，通常也使用后者更加合理1.2误差基础知识十进制数的标准形式(其中),1.2.3有效数字四舍五入保留位:因此有误差限称为有效数字,

称为有效数问题:有效数字和相对误差限有什么关系?有效数的误差限是末位数单位的一半,其本身就体现了误差界有效数末尾不可以随便增减零1.2误差基础知识1.2.4计算机的浮点数系单精度实数由32位二进制的浮点数表示:1位符号,

23位尾数,

8位阶数(本身也有符号)机器所能表示的数中,离最近的是和.因此若,则在机器中记为,即.相对误差限最大数,最小数,

上溢,

下溢1.2误差基础知识1.2.5一个实例有一艘驳船,宽度为5米,欲驶过一个河渠.该河渠有一个直角弯道,形状和尺寸如图所示.试问,要驶过这个河渠,驳船的长度不能超过多少米?驳船的长度有如下关系极小化问题1.2误差基础知识1.2.5一个实例这个过程中有多少处有误差？或者求零点问题可证,对任意并可求得1.2误差基础知识1.2.6数值计算中应注意的几个问题计算容易推导出因此,调用MatLab程序nademo11.2误差基础知识1.2.6数值计算中应注意的几个问题计算,其中,保留四位有效数字.其它的例子:避免大数和小数相加减1.2误差基础知识1.2.6数值计算中应注意的几个问题,其中计算中间步骤简化计算步骤Horner算法或秦九韶算法1.3Matlab软件1.3.1简介全称:MatrixLaboratory功能:科学计算、符号计算、图形处理等数据类型:数、字符串、矩阵、单元型数据和结构型数据集成界面:命令窗口、命令历史窗口、当前路径窗口、工作空间变量窗口等提示符>>,换行符...,注释符

%,默认变量ans1.3Matlab软件1.3.2向量和矩阵的基本运算矩阵A=[13;24]向量a=[123456]冒号a=1:6a:s:b列向量A=[1;2;3]字符串A='hellomatlab'

A='This''smatlab''sworld.'1.3Matlab软件1.3.2向量和矩阵的基本运算常量:在运行过程中不能变化的量科学记数法:显示方式:format(只影响显示)变量:保存在内存(地址),可随时变化内置变量:i,j,pi,Inf,NaN(NotaNumber)Inf及NaN的运算规律1.3Matlab软件1.3.2向量和矩阵的基本运算矩阵的加(+)、乘(*)、数乘(*)、幂(^)矩阵的点乘(.*)、点除(./)、点幂(.^)矩阵的左除(X=A\B即求解AX=B)矩阵的右除(X=A/B即求解A=XB)1.3Matlab软件1.3.2向量和矩阵的基本运算>>x=[0pi/6pi/4pi/3pi/2];>>sin(x)ans=

0

0.5000

0.7071

0.8660

1.0000向量功能其他初等函数:三角反三角、指数对数、根号、绝对值等等1.3Matlab软件1.3.2向量和矩阵的基本运算>>sqrt([91011])>=pians=

0

1

1>>a=[2300];>>b=[-1010];>>a&bans=

1

0

0

0>>a|bans=

1

1

1

0>>~bans=

0

1

0

11.3Matlab软件1.3.2向量和矩阵的基本运算矩阵运算>>A=magic(3)

816>>A(2,1:3)

357>>A(2:end,[1end])492

>>B=[23];>>C=[12;34];>>D=[57]';>>A=[B9;CD]

>>A(A>=4)=0

>>v=1:9;>>v(abs(v-5)<=2)=[]1.3Matlab软件1.3.3流程控制基本语法ifvalue1,

statement1,elseifvalue2,

statement2,else

statement3end1.3Matlab软件1.3.3流程控制例如(判别闰年)ifmod(year,400)==0,

fprintf('%disaleapyear.\n',year);elseifmod(year,100)==0,

fprintf('%disnotaleapyear.\n',year);elseifmod(year,4)==0,

fprintf('%disaleapyear.\n',year);else

fprintf('%disnotaleapyear.\n',year);end1.3Matlab软件1.3.3流程控制基本语法forloopvalue=value,

statement,end和whilevalue,

statement,end1.3Matlab软件1.3.3流程控制例如(利用计算圆周率的近似值)>>s=0;>>fork=1:10000,

s=s+1/k^2;

end>>s=sqrt(6*s)1.3Matlab软件1.3.3流程控制Collatz猜想：输入一个正整数n,如果是偶数就除以2,是奇数就乘3加1,如此一直变换,最后会变成1.n=input('n=

');whilen~=1,

ifmod(n,2)==1,

n=n*3+1;

else

n=n/2;

end

disp(n);end1.3Matlab软件1.3.3流程控制冒泡排序：把一列数想象为垂直存放,数值大的在下方,每轮比较时从上到下依次比较相邻的两个数,若是上面的数大,把它们对调,否则不动。直至没有对调为止。>>done=0;k=1;>>v=input('arowvector:');arowvector:

[1863975024]>>while~done,

done=1;

forp=1:length(v)-k,

ifv(p)>v(p+1),

tmp

=v(p);

v(p)=v(p+1);

v(p+1)=tmp;

%OR

v([pp+1])=v([p+1p]);

done=0;

end

end

k=k+1;

end1.3Matlab软件1.3.4脚本文件和函数文件把一系列命令收集在一个文件里,保存为以.m为后缀的文件.执行时只需要键入文件名,不需键入后缀.例:>>mysortarowvector:

[1863975024]v=

0

1

2

3

4

5

6

7

8

91.3Matlab软件1.3.4脚本文件和函数文件一种封装的文件,具有特定的头格式:function[out1,out2,...]=funname(in1,in2,...)函数名必须和文件名一致与脚本文件的比较例:文件mysort2.m1.3Matlab软件1.3.4脚本文件和函数文件函数头function[v,s]=mysort3(v)调用>>d=[53421];>>[r,w]=mysort3(d)传值方式:在输出或输入列表中的位置列表不一样长的情形命令global的用法1.3Matlab软件1.3.4脚本文件和函数文件变量nargin和nargout的含义用法(例如:根据输入计算面积)functions=zhouchang(a,b,c)ifnargin==1,

s=2*pi*a;elseifnargin==2,

s=2*(a+b);elseifnargin==3,

s=a+b+c;end1.3Matlab软件1.3.4脚本文件和函数文件直接或间接地用到了自己例如:Fibonacci数列定义为functionf=fib(n)ifn>=3,

f=fib(n-1)+fib(n-2);elseifn==1|n==2,

f=1;end1.3Matlab软件1.3.5帮助系统Help：查看工具箱,函数可以自己书写文件的帮助,写在function之后其他查看系统命令用法的工具:doc,lookfor其他帮助命令:which,who等辅助命令:clc,home,clear1.3Matlab软件1.3.6画图功能>>x=0:0.01:10;>>y=1./(1+x.^2)+sin(x).*exp(x/3);plot(x,y,‘g*-’)画函数的图像hold命令命令plot中的参数选项plot(x,y1,'yo--',x,y2,'g*',x,y3,'r+:',x,y4,'bp:’);其他类似命令:loglog,semilogx,semilogy1.3Matlab软件1.3.6画图功能三维线图>>t=linspace(0,10*pi,2000);>>plot3(sin(t).*t,cos(t).*t,t,'r-');>>view(-17,66)1.3Matlab软件1.3.6画图功能三维面图：命令meshgrid>>x=1:4;>>y=5:3:11;>>[X,Y]=meshgrid(x,y)X=

1

2

3

4

1

2

3

4

1

2

3

4Y=

5

5

5

5

8

8

8

8

11

11

11

111.3Matlab软件1.3.6画图功能三维面图：命令surf及contour例如:画下面函数的图像及等高线>>x=linspace(-10,10,200);>>[X,Y]=meshgrid(x);>>Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);>>surf(X,Y,Z);>>contour(X,Y,Z,20)1.3Matlab软件1.3.6画图功能标注:坐标轴xlabel,曲线legend,图形标题title窗口控制:打开figure,关闭close,清除clf坐标轴控制:axis('equal','off',[-1327])或者

axisequal;

axisoff1.3Matlab软件1.3.7数据操作文本方式>>x=0:0.1:1;y=[x;exp(x)];>>fid=fopen('exp.txt','wt');>>fprintf(fid,'%s\n','%---exp.txt---');>>fprintf(fid,'%6.2f%12.8f\n',y);>>fclose(fid);>>fid=fopen('exp.txt');>>s=fscanf(fid,'%c',[117])s=%---exp.txt--->>A=fscanf(fid,'%6f%12f\n',[2inf]);>>A=A'转置关系1.3Matlab软件1.3.7数据操作二进制方式>>loadclown.mat>>whoYourvariablesare:X

caption

map>>image(X);colormap(map)>>savea.matx*A1.3Matlab软件1.3.7数据操作>>A=[120;030;0-16];>>A=sparse(A)A=

(1,1)

1

(1,2)

2

(2,2)

3

(3,2)

-1

(3,3)

6或者1.3Matlab软件1.3.7数据操作>>i=[11233];>>j=[12223];>>s=[123-16];>>A=sparse(i,j,s)A=

(1,1)

1

(1,2)

2

(2,2)

3

(3,2)

-1

(3,3)

6指明阶数>>A=sparse(i,j,s,200,100);找出非零元>>[i,j,s]=find(A);1.3Matlab软件1.4评注Matlab参考书目:[1]MatLab与科学计算(第2版),王沫然,电子工业出版社,2007年8月[2]MatLab工程数学应用,许波、刘征,清华大学出版社,2000年4月[3]MatLab数学实验,胡良剑、孙晓君,高等教育出版社,2006年6月[4]MatLab高等数学实验,章恩栋、马玉兰、徐美萍、李双,电子工业出版社,2008年11月AdvancedNumericalComputing现代数值计算同济大学数学科学学院学海无涯，祝你成功！第二章AdvancedNumericalComputing线性代数方程组的直接法现代数值计算同济大学数学科学学院目录/Contents第二章线性代数方程组的直接法第一节高斯消去法第二节矩阵的三角分解第三节正交矩阵与奇异值分解前言求解大规模线性方程组:如何利用计算机来快速、稳定、有效地求解该问题是科学计算的核心问题之一直接法和迭代法由于浮点运算的精度的影响，直接法不可能给出完全精确的计算解给定

阶方阵

和

维向量

,寻找向量

使得：2.1高斯消去法2.1.1基本步骤记线性方程组的分量形式为下面演示用高斯消去法求解上述线性方程组的计算过程。计算过程分为消去过程和回代过程记矩阵

，向量

，它们的元素分别为：一、消去过程第一步：如果

，用数

依次乘以方程组的第一行，并加到第

行上去，

其中2.1高斯消去法2.1.1基本步骤第二步：如果

，用数

依次乘以方程组的第二行，并加到第

行上去，其中2.1高斯消去法2.1.1基本步骤类似地，这样的运算过程一直可作到第

步，结果转化为一个上三角形方程组2.1高斯消去法2.1.1基本步骤二、回代过程如果

，可以逐次回代计算出线性方程组的解。这就是求解线性方程组的高斯消去法。在没有浮点运算误差的情况下，该方法在有限的计算步骤内能够得到原线性方程组的精确解，是直接法的一种。2.1高斯消去法2.1.1基本步骤【高斯消去法】算法(1)对

做：

对

做：

用数

乘以方程组第

行加到第

行上;

标记得到的矩阵及右端向量的上标 ;(2) 且对于

做:2.1高斯消去法2.1.1基本步骤消去过程的第

步，对矩阵需作

次乘法运算及

次除法运算，对右端向量作

次乘法运算，在消去过程总的乘除法运算工作量为回代过程中，计算每个

需作

次乘除法运算，其工作量为用高斯消去法计算线性方程组所需要总的乘除法运算工作量为2.1高斯消去法2.1.2乘除法的运算量高斯消去法能够顺利进行到底是有前提条件的，即要求所有的主导元素不等于零。如果某个主元为零，则高斯消去法中断。【定义2.1】在计算中做除数的元素

，被称为主导元素，简称主元。2.1高斯消去法2.1.3选主元策略【例2.1】取

，用高斯消去法计算下述线性方程组。（假定模型计算机具有8位字长的浮点表示及16位的累加器。）解：首先用高斯消去法，这时

，对方程组消元，于是得到2.1高斯消去法2.1.3选主元策略在这个模型计算机上，具体计算是这样的：

通过四舍五入，输出结果变成

同理，

的计算结果也是

。这样高斯消去的计算解为：

实际上方程组的精确解为

这里未知量

的计算解的相对误差达到了惊人的

，这就是“大数吃小

数”的现象。2.1高斯消去法2.1.3选主元策略由于浮点运算误差的影响，高斯消去法过程中会得到错误的解。为了避免上述不稳定的现象，对一般的线性方程组而言，我们采用选主元的策略。采用列主元素高斯消去法，对上述例子中方程组进行行交换：这时

，进行消元后，可得到在模型计算机上，

和

都被算成为

。2.1高斯消去法2.1.3选主元策略若考虑方程此时,选不选主元一样:因为这种方程的选主元同等变形变得毫无意义。两边同乘了一个很大的数2.1高斯消去法2.1.3选主元策略目录/Contents第一节高斯消去法第二节矩阵的三角分解第三节正交矩阵与奇异值分解第二章线性代数方程组的直接法2.2矩阵的三角分解常见的矩阵三角分解有：LU分解（杜利脱尔分解，克洛脱分解）LDU分解乔列斯基分解【定义2.1】把一个

阶矩阵分解成结构简单的三角形矩阵的乘积称为矩阵的三角分解。高斯消去法的消去过程与左乘下述矩阵是等价的：2.2.1LU分解2.2矩阵的三角分解其中易知则其中2.2.1LU分解2.2矩阵的三角分解这里，是单位下三角矩阵，

是上三角矩阵，这种矩阵分解称为杜利脱尔(Doolittle)分解，或者杜利脱尔三角分解。克洛脱(Crout)分解LDU分解以上三种分解统称为矩阵的三角分解，或者

分解。如果不作特殊说明，一般我们所说的

分解就是指杜利脱尔三角分解。，这里是下三角矩阵，

是单位上三角矩阵。

，这里

是单位下三角矩阵，

是对角矩阵，是单位上三角矩阵。2.2.1LU分解2.2矩阵的三角分解矩阵三角分解的存在唯一性【定理2.1】(存在性)利用高斯消去法求解方程组

时的主元素

的充要条件是

阶矩阵

A

的所有顺序主子式均不为零。若

A

为

阶矩阵，且所有顺序主子式均不等于零，则

A

可分解为一个单位下三角矩阵

L与一个上三角矩阵

U

的乘积: ，且分解是唯一的。【定理2.2】(唯一性)2.2.1LU分解2.2矩阵的三角分解【算法2.1】杜利脱尔算法(1)对

做:(2)(3)2.2.1LU分解2.2矩阵的三角分解(1)对

做:(2)(3)2.2.1LU分解2.2矩阵的三角分解【算法2.2】克洛脱算法【算法2.3】回代(1)对

做：(2)对

做：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩阵的三角分解解：【例2.1】利用LU分解求解

：用MATLAB可以计算矩阵的LU分解，其语法为:>>A=[211;

232;

234];

>>[l,u]=lu(A)

l=

1

0

0

1

1

0

1

1

1

u=

2

1

1

0

2

1

0

0

22.2.1LU分解2.2矩阵的三角分解>>b=[479]';>>y=L\b;

y=

4

3

2>>x=U\y

x=

1

1

12.2.1LU分解2.2矩阵的三角分解【例2.2】已知

，作A的杜利脱尔分解，并求解方程组，其中解：假设2.2.1LU分解2.2矩阵的三角分解经计算可得对

进行回代可得

；再对

进行回代，则

。用MATLAB可以计算矩阵的LU分解，其语法为:2.2.1LU分解2.2矩阵的三角分解[L,U]=lu(A)当矩阵A为对称正定时，它的所有顺序主子式都大于零，易知存在唯一的LU分解：由A的对称性可得

按照分解的唯一性可得

即

2.2.2乔列斯基分解2.2矩阵的三角分解易证这个三角分解是唯一的，称之为乔列斯基(Cholesky)分解根据

A

是对称正定矩阵，有从而

记

则有其中，

G是对角元均大于零的下三角矩阵2.2.2乔列斯基分解2.2矩阵的三角分解【算法2.4】乔列斯基分解算法(1)对

(2) (3) 2.2.2乔列斯基分解2.2矩阵的三角分解给定乔列斯基分解，线性方程组

的求解可转化为计算公式为平方根法2.2.2乔列斯基分解2.2矩阵的三角分解解：设A的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组2.2.2乔列斯基分解2.2矩阵的三角分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

A

的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组2.2.2乔列斯基分解2.2矩阵的三角分解解：设

A

的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组2.2.2乔列斯基分解2.2矩阵的三角分解解：设

A的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组2.2.2乔列斯基分解2.2矩阵的三角分解解：设A的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组2.2.2乔列斯基分解2.2矩阵的三角分解解：设

A

的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组回代2.2.2乔列斯基分解2.2矩阵的三角分解解：设

A

的乔列斯基分解

，经计算得【例2.3】利用平方根法求解下述对称正定方程组【例2.4】求下列矩阵的乔列斯基分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得【例2.4】求下列矩阵的乔列斯基分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得【例2.4】求下列矩阵的乔列斯基分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得【例2.4】求下列矩阵的乔列斯基分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得【例2.4】求下列矩阵的乔列斯基分解2.2.2乔列斯基分解2.2矩阵的三角分解解：设

的乔列斯基分解

，经计算得【例2.4】求下列矩阵的乔列斯基分解利用矩阵的三角分解，很容易导出一些特殊方程组的解法。对矩阵

作克洛脱分解，得2.2.3追赶法2.2矩阵的三角分解设

为三对角矩阵，即(3）(4)(1)设

(2)对

算法2.52.2.3追赶法2.2矩阵的三角分解追赶法的Matlab程序tridiagsolver.m如下functionx=tridiagsolver(A,b)[n,n]=size(A);fori=1:nif(i==1) l(i)=a(i,i); y(i)=b(i)/l(i); else l(i)=a(i,i)-a(i,i-1)*u(i-1); y(i)=(b(i)-y(i-1)*a(i,i-1))/l(i); end if(i<n) u(i)=a(i,i+1)/l(i); end2.2.3追赶法2.2矩阵的三角分解endx(n)=y(n)forj=n-1:-1:1 x(j)=y(j)-u(j)*x(j+1);end2.2.3追赶法2.2矩阵的三角分解追赶法的Matlab程序tridiagsolver.m如下【例2.5】用追赶法求解下述三对角线性方程组解：追的过程为2.2.3追赶法2.2矩阵的三角分解赶的过程为因此，原方程组的解为

。2.2.3追赶法2.2矩阵的三角分解设满足且

和

分别是分块下三角矩阵和分块下三角2.2.3分块三角分解2.2矩阵的三角分解矩阵经计算即这样的矩阵分解称为矩阵的分块三角分解2.2.3分块三角分解2.2矩阵的三角分解其中，

称为

的舒尔(Schur)补【例2.6】求分块矩阵的一个分块三角分解，其中解：因为所以，2.2.3分块三角分解2.2矩阵的三角分解因此，

分块矩阵

的一个分块三角分解为:2.2.3分块三角分解2.2矩阵的三角分解目录/Contents第一节高斯消去法第二节矩阵的三角分解第三节正交矩阵与奇异值分解第二章线性代数方程组的直接法矩阵除了三角分解以外，还有QR分解为正交矩阵，

为上三角矩阵。奇异值分解和

为正交矩阵，

为对角矩阵。2.3正交矩阵与奇异值分解2.3.1正交矩阵正交矩阵

有如下性质：

的长度与

的长度相等矩阵

，满足

，

我们称这样的矩阵

为正交矩阵。2.3正交矩阵与奇异值分解单位矩阵任意个置换矩阵的乘积仍然是正交矩阵置换矩阵：将单位矩阵的任意两行（列）交换得到的矩阵。

譬如，交换第

行和第

行，得：2.3.1正交矩阵2.3正交矩阵与奇异值分解旋转矩阵（Givens变换）形如的矩阵被称为Givens矩阵或Givens变换，或称（平面）旋转矩阵（或旋转变换）其中

为旋转的角度。2.3.1正交矩阵2.3正交矩阵与奇异值分解平面旋转变换的几何意义：2.3.1正交矩阵2.3正交矩阵与奇异值分解反射矩阵（Householder变换）设

，且

，则称为Householder变换，或者Householder矩阵。

2.3.1正交矩阵2.3正交矩阵与奇异值分解

Householder矩阵有如下性质：不难验证所以，Householder变换又称镜面反射变换，Householder矩阵也称初等反射矩阵。(1)

，即

是对称阵。(2)

，即

是正交阵。(3)设

，

为过原点的平面且

，可分解成为

，其中

。2.3.1正交矩阵2.3正交矩阵与奇异值分解镜面反射变换的几何意义：2.3.1正交矩阵2.3正交矩阵与奇异值分解则从而其中,其中，

。一个重要的应用是对

，求Householder矩阵

，使得由

知

，由

的构造，有设

，为了使

计算时不损失有效数位，取2.3.1正交矩阵2.3正交矩阵与奇异值分解则【例3.1】已知

，求Householder矩阵

，使得

，

其中

解：取

，

，

，2.3.1正交矩阵2.3正交矩阵与奇异值分解证：给出构造性证明如下：【定理3.1】设

，则存在正交阵

，使

，其中

为上三角阵。首先，考虑

的第一列

，可找到Householder矩阵

，使得

的除了第

个元素以外都为零。

同理，找到

使得

的第2列对角元以下元素为0，而第一列对角元以下元素与

一样是0。依次这样下去，可以得到其中

为上三角矩阵，

为正交阵。定理证毕。2.3.2QR分解2.3正交矩阵与奇异值分解【定理3.2】

，设

非奇异，则存在正交阵

与上三角阵

，使得

有如下分解：且当

的对角元均为正时，分解是唯一的。该定理保证了

可分解为

，若

非奇异，则

也非奇异。如果不规定

的对角元为正，则分解不是唯一的。2.3.2QR分解2.3正交矩阵与奇异值分解【例3.2】用Householder变换作矩阵

的QR分解则有解：Householder矩阵

，使2.3.2QR分解2.3正交矩阵与奇异值分解和再找

，使

，得且这是一个下三角矩阵，但对角元皆为负数。2.3.2QR分解2.3正交矩阵与奇异值分解只要令

，则有

是对角元为正的上三角矩阵，使得

。其中，2.3.2QR分解2.3正交矩阵与奇异值分解

【定理3.3】（奇异值分解）设

，则存在酉阵

和

，使得：其中

, , 。【定义3.2】设

，

的

个特征值的非负平方根叫作

的奇异值，记为

。2.3.3奇异值分解2.3正交矩阵与奇异值分解奇异值分解与对称矩阵基于特征向量的对角化类似，但还是有明显的不同。对称矩阵特征向量分解的基础是谱分析，而奇异值分解则是谱分析理论在任意矩阵上的推广。奇异值分解提供了一些关于

的信息，例如非零奇异值的数目（

的阶数）和

的秩相同。一旦秩r确定，那么U的前r列构成了

的列向量空间的正交基。奇异值分解非常有用和可靠的分解，但是它比QR分解要花上近十倍的计算时间。2.3.3奇异值分解2.3正交矩阵与奇异值分解本章小结高斯消去法矩阵的三角分解法矩阵的QR分解法直接法相对来说，工作量小，精度高，但程序复杂，并且对于高阶线性方程组易于受计算机容量的限制，所以它适于求解中小型方程组。对于高阶大型线性方程组，有效的解法是第三章要讨论的迭代法。AdvancedNumericalComputing现代数值计算同济大学数学科学学院学海无涯，祝你成功！第三章AdvancedNumericalComputing线性代数方程组的迭代法现代数值计算同济大学数学科学学院目录/Contents第三章线性方程组的迭代解法第一节范数和条件数第二节基本迭代法第三节不定常迭代法迭代方法的优点：能够充分利用系数矩阵的稀疏性质占用内存少运算方便，计算程序也较为简单前言迭代法分类常用定常迭代法雅可比（Jacobi）迭代法高斯-塞德尔(Gauss-Seidel)迭代法超松弛(SOR)迭代法常用不定常迭代法共轭梯度法广义极小残量（GMRES）法前言前言基本概念【定义3.1】线性方程组的解是一个向量，称为解向量【定义3.2】近似解向量与精确解向量之差称为近似解的误差向量为了估计误差的大小，需要引入衡量向量和矩阵大小的度量概念-范数，满足：【定义3.4】对任意

阶方阵

，若对应一个非负实数时成立；（1）

，等号当且仅当（2）对任意实数，

；（3）对任意两个n阶方阵

和

，

；（4）对任意两个n阶方阵

和

，

；则称

为方阵

的矩阵范数。为的谱半径3.1范数和条件数3.1.1向量范数和矩阵范数【定义3.5】设

是

阶矩阵，则称这里，

为

的特征值（

）对于

阶方阵，定义分别为矩阵的1范数，2范数，无穷范数和范数。3.1范数和条件数3.1.1向量范数和矩阵范数【定义3.6】满足的向量范数和矩阵范数，称为相容的。常用范数有下列相容关系：3.1范数和条件数3.1.1向量范数和矩阵范数考虑线性方程组解是由

和

决定的；数据

和

都带有误差（测量误差，模型误差，舍入误差）；3.1范数和条件数3.1.2扰动分析和条件数通常

和

的误差相对于精确数据都是微小的，称之为小的扰动，分别记为

和

。一个很自然的问题是：

和

的微小扰动将对计算线性方程组的解有何影响？考察解为若

，则原方程组变为解为3.1范数和条件数3.1.2扰动分析和条件数有即解的相对误差是右端项相对误差的10000倍。表明一些线性方程组系数的微小变化会引起解的巨大变化。为什么呢？3.1范数和条件数3.1.2扰动分析和条件数若

有一个小扰动

，则解

产生一个扰动

，即：

有

两边取范数另得：3.1范数和条件数3.1.2扰动分析和条件数若

有一个小扰动

，则解

相应产生一个扰动

，即：有两边取范数得得3.1范数和条件数3.1.2扰动分析和条件数观察得，解

的相对误差除了与方程组系数矩阵

和右端

扰动有关外，还与

的大小有关。条件数的大小与范数有关：【定义3.7】设

为

阶非奇异矩阵，称数

为线性方程组

的条件数，或者称为矩阵

的条件数。3.1范数和条件数3.1.3向量范数与矩阵范数如果

对称正定，则和

是最大和最小特征值。条件数有如下性质：3.1范数和条件数3.1.3向量范数与矩阵范数对于任意的

阶非奇异矩阵

，

成立；对于任意的正交矩阵

，有

；对于任意的

阶非奇异矩阵

及任意非零常数

，

成立；对于任意的

阶非奇异矩阵

及任意

阶正交矩阵

，有典型的病态矩阵希尔伯特（Hilbert）矩阵随着

的增大,条件数

非常迅速的增加希尔伯特矩阵阶数越高，病态程度就越为严重。3.1范数和条件数3.1.3向量范数与矩阵范数目录/Contents第三章线性方程组的迭代解法第一节范数和条件数第二节基本迭代法第三节不定常迭代法假定

有分解为

；其中

是非奇异方阵3.2

基本迭代法考虑有或；其中

且

3.2.1基本知识【算法】建立迭代公式：即给定初始向量

，按上述公式迭代，可以得到一个

向量序列

这种方法就是解线性方程组的基本迭代解法若

收敛于确定的向量

，则

，亦即

，

就是方程组

的解3.2

基本迭代法3.2.1基本知识令其中为对角矩阵，严格下三角矩阵和严格上三角矩阵通过不同的构造方法，可得到三种基本迭代解法：雅可比迭代法高斯-塞德尔迭代法SOR迭代法3.2

基本迭代法3.2.1基本知识线性代数方程组Jacobi迭代:从第

个方程解出

3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法Jacobi迭代:从第

个方程解出

即取

和

，得雅可比迭代法:或：迭代矩阵为：【算法3.1】雅可比迭代算法1)选定初值2)对

计算:3)如果近似解达到收敛条件，退出；否则，继续第2)步的计算。3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法function[x,iter]=jacobi(A,b,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=D\(b+L*x+U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法Gauss-Seidel迭代:从第

个方程解出 ,尽量使用新分量3.2

基本迭代法3.2.3高斯-塞德尔(Gauss-Seidel)迭代法如果雅可比迭代法是收敛的，将

的分量投入下一个迭代方程中使用，这样可能会收到更好的效果，即高斯-塞德尔迭代法:或迭代矩阵为：3.2

基本迭代法3.2.3高斯-塞德尔(Gauss-Seidel)迭代法【算法3.2】高斯-塞德尔迭代算法3)如果近似解达到收敛条件，退出；否则，继续第2)步的计算1)选定初值

2)对

计算:3.2

基本迭代法3.2.3高斯-塞德尔(Gauss-Seidel)迭代法function[x,iter]=gs(A,b,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=(D-L)\(b+U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.3高斯-塞德尔(Gauss-Seidel)迭代法迭代格式为:即：为了加快迭代的收敛速度，等号右端的第二项前乘以一个参数

，得逐次超松弛迭代法，简称SOR迭代法，

被称为松弛因子3.2

基本迭代法3.2.4超松弛(SOR)迭代法迭代矩阵为：注意当

时，SOR迭代法其实就是GS迭代法【算法3.3】SOR迭代算法1)选定初值

,3)如果近似解达到收敛条件，退出；否则，继续第2)步的计算。2)对

计算3.2

基本迭代法3.2.4超松弛(SOR)迭代法function[x,iter]=sor(A,b,omega,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=(D-omega*L)\(b+(1-omega)*D*x+omega*U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.4超松弛(SOR)迭代法【例3.1】取初值为

，试用雅可比迭代，GS迭代以及 SOR迭代分别计算线性方程组迭代过程保留5位有效数字（精确解

）3.2

基本迭代法3.2.4超松弛(SOR)迭代法建立雅可比迭代格式解：雅可比方法迭代21步就得到了保留5位有效数字的近似解取

，

对

计算可得

3.2

基本迭代法3.2.4超松弛(SOR)迭代法建立GS迭代格式取

，对

计算可得GS方法迭代9步就得到了保留5位有效数字的近似解。对于这个例子而言，GS迭代的收敛速度大约是雅可比迭代的2倍。3.2

基本迭代法3.2.4超松弛(SOR)迭代法选取

，建立SOR迭代格式:取

，对

计算可得3.2

基本迭代法3.2.4超松弛(SOR)迭代法SOR方法迭代7步就得到了保留5位有效数字的近似解。对于这个例子而言，

时，SOR迭代法的收敛速度比GS迭代法要快。首先引进两个重要的概念:【定义3.8】设

是

阶矩阵

，如果存在

阶排列矩阵

，

使

，有如下形状：3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计其中，

和

分别为

阶和

阶方阵

，则称

为可约矩阵，如果不存在这样的排列阵，则称

为不可约阵。首先引进两个重要的概念:【定义3.9】设

是

阶矩阵，若

满足：且其中至少有一个严格不等式成立，称

是(行)弱对角占优若式中的每一个不等式都是严格不等号，则称

是(行)严格对角占优3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计【引理3.1】设

阶矩阵

是严格对角占优或不可约弱对角占优矩阵，

则

是非奇异矩阵。证明(一)：假设

不可逆。则存在非零向量

，使得

不妨设

且

则

从

的第

个方程得是严格占优的，因此上述不等式不成立，矛盾。3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计利用矩阵的Jordan标准型，可以得出下述引理【定理3.10】任一矩阵

的谱半径均不大于

的任一与某一向量范数相容的

矩阵范数，即证明:两边取范数，假设

是矩阵

的特征对，则【引理3.2】设

是任意

阶矩阵，则

的

次幂

（当

）的

充要条件为谱半径3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计【定理3.11】对于基本迭代格式，给定初值

，

是真解，有下列收敛

结果和误差估计：(1)迭代格式收敛的充要条件为谱半径

；(2)如果

，则有如下估计：3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计证明:是任意的，上式等价于

，即 .3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计【定理3.11】对于基本迭代格式，给定初值

，

是真解，有下列收敛

结果和误差估计：(1)迭代格式收敛的充要条件为谱半径

；证明:移项即得3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计【定理3.11】对于基本迭代格式，给定初值

，

是真解，有下列收敛

结果和误差估计：(2)如果

，则有如下估计：【定理3.12】若

是严格对角占优或不可约弱对角占优矩阵，则雅可比迭代和 GS迭代都收敛【定理3.14】SOR迭代收敛的必要条件是

【定理3.13】若

是对称正定矩阵，则雅可比迭代收敛的充要条件是

也

对称正定【定理3.15】设系数矩阵

对称正定，则

时SOR迭代收敛3.2

基本迭代法3.2.5迭代的收敛性分析和误差估计目录/Contents第三章线性方程组的迭代解法第一节范数和条件数第二节基本迭代法第三节不定常迭代法本节将介绍两类最基本的不定常迭代方法：一类是求解对称正定线性方程组的最速下降法和共轭梯度法另一类是求解不对称线性方程组的广义极小残量法3.3

不定常迭代法对于一切

，通过直接计算：知函数

的梯度为：定义

元二次函数

为：3.3

不定常迭代法3.3.1最速下降法令

为方程组的真解，则有对于一切

，

，且对于一切

有

3.3

不定常迭代法3.3.1最速下降法【定理3.16】设

对称正定，

是方程组

的解的充要条件是

为

二次函数的极小值点，即：

设

，则由

的对称正定性有：上式对于任意的

都成立，即

充分性设

，则对任意的

，可知

，对任意的向量

都成立。由此可得

成立3.3

不定常迭代法3.3.1最速下降法选定初值

，对

直到收敛，计算

结束算法【算法3.4】最速下降法3.3

不定常迭代法3.3.1最速下降法【引理3.3】设

是

阶实对称正定矩阵，

是

次实系数多项式，

则对任意

有【定理3.17】设

是

阶实对称正定矩阵，

和

分别是

的最大和最小

的特征值，则由最速下降法得到的迭代序列

满足下述

误差估计：3.3

不定常迭代法3.3.1最速下降法

其中，

是

的特征值，

是

中的向量范数对上述最速下降法作一简单分析，就会发现，负梯度方向虽从局部来看是最佳的搜索方向，但从整体来看并非最优。共轭梯度法还是采用一维极小搜索的思想，不再沿负梯度方向方向搜索，而要另找一组新的方向进行

次一维搜索后，求得近似解【定义3.18】

对称正定，若

中向量组

满足则，称它为

中的一个A-共轭向量组，或者称A-正交向量组3.3

不定常迭代法3.3.2共轭梯度法【算法3.5】共轭梯度法选定初值

，设

，对

直到收敛，计算

如果

足够小退出算法；

否则

结束算法3.3

不定常迭代法3.3.2共轭梯度法下面不加证明地给出一些理论结果：【定理3.19】由共轭梯度法得到的向量组

和

具有如下性质:1)2)3)【定理3.20】设

是

阶实对称正定矩阵，

和

分别表示

的最大和

最小特征值，则由共轭梯度法得到的向量序列

满足下述

误差估计：3.3

不定常迭代法3.3.2共轭梯度法【例3.2】试用共轭梯度法求解下述线性方程组，其中，初值为

，使得最终迭代误差

达到

，最大迭代步数设为10。3.3

不定常迭代法3.3.2共轭梯度法设，则：解：第一步迭代如下：我们发现

3.3

不定常迭代法3.3.2共轭梯度法因此进行第二步迭代如下：此时

，满足迭代收敛条件，因此迭代终止。综上，总共需要两步迭代达到其中迭代解为

。3.3

不定常迭代法3.3.2共轭梯度法

广义极小残量法

（GerneralMinimalRESidual)

这个算法经二十世纪八十年代来许多研究者的努力得到很大的完善，已经成为当前求解大型稀疏非对称线性方程组的主要手段。我们从

开始，构造一组相互正交且范数为1的向量

,使得其中为一个上海瑟伯格(Hessenberg)矩阵3.3

不定常迭代法3.3.3广义极小残量法令有假设存在一个

满足

，我们可以看到：其中，

是第一个分量为1，其余分量为0的

维列向量3.3

不定常迭代法3.3.3广义极小残量法选定初值

，设

和

对

直到收敛，计算

对

，计算

【算法3.6】广义极小残量法3.3

不定常迭代法3.3.3广义极小残量法如果

足够小退出算法；

否则

求解最小二乘问题

得到的解为

计算

结束算法下面不加证明地给出广义极小残量法的收敛性定理：假设

是可对角化的，即

，

，且

的全部特征值落在不包含原点的椭圆

内部，

分别代表椭圆的中心，焦距和长半轴。则，

其中，

是

阶的Chebyshev多项式【定理3.21】3.3

不定常迭代法3.3.3广义极小残量法预处理技术从广义上来说可以指对原方程组进行的任何显示的或者隐式的修正，使得该方程组通过迭代法更容易求解。分别称为左预处理或者右预处理的迭代方法。譬如，找到一个矩阵

，求解3.3

不定常迭代法3.3.4预处理技术这里，

被称为预处理矩阵。一个好的预处理矩阵至少能够满足如下的条件：构造

的代价很小；跟

足够接近；关于

的线性方程组很容易求解。3.3

不定常迭代法3.3.4预处理技术

总结本节介绍了一些常用的迭代法及其收敛性质:雅可比（Jacobi）迭代法高斯-塞德尔(Gauss-Seidel)迭代法超松弛(SOR)迭代法共轭梯度法

广义极小残量法AdvancedNumericalComputing现代数值计算同济大学数学科学学院学海无涯，祝你成功！第四章AdvancedNumericalComputing多项式插值与样条插值现代数值计算同济大学数学科学学院目录/Contents第四章多项式插值与样条插值第一节多项式插值第二节拉格朗日（Lagange）插值第三节牛顿（Newton）插值第四节埃尔米特（Hermite）插值第五节三次样条插值前言本章和下一章论述的主题：函数表达的问题第一个问题，假定已有一个函数的数值表(表1).提问：是否能找到一个简单而便于计算的公式，利用它可以精确地重新算得这些给定点。------插值问题表1:

......第二个问题，设给定一个函数

，

的表达式非常复杂，计算

的值很不经济。在这种情况下，就要寻找另一个函数

，它既易于求值且又是对

的一个合理的逼近。——连续函数的逼近问题前言第三个问题，假定表中给出的数值带有误差。比如当这些值来自于物理实验时，就可能出现这种情况。现在要寻找一个公式，使得它可以近似地表示这些数据。

——离散函数的逼近问题前言4.1多项式插值4.1.1多项式插值问题的定义【定义4.1】设函数

在区间

上有定义，且已知它在

个互异的点

的函数值

若存在一个次数不超

次的多项式

(其中

为实数)满足条件：则称

为函数

的

次插值多项式按上述条件求函数

的近似表达式

的方法称为插值多项式插值法，称该条件为插值多项式插值条件, 为插值多项式插值节点，包含插值节点的区间

为插值多项式插值区间几个常用术语:插值法------按插值条件求函数

的近似表达式

的方法插值条件------条件(2.1)插值节点------插值区间------包含插值节点的区间

4.1多项式插值4.1.1多项式插值问题的定义多项式插值的几何意义：4.1多项式插值4.1.1多项式插值问题的定义4.1多项式插值4.1.1多项式插值问题的定义多项式插值的几何意义：4.1多项式插值4.1.1多项式插值问题的定义多项式插值的几何意义：证：设

次多项式

是函数

在

上的

个互异的节点

的插值多项式，则求

的问题就

可归结为求它的系数

为未知元的

阶线性方程组：其系数行列式是范德蒙德（Vandermonde)行列式：【定理4.1】

次插值多项式存在且唯一（4.1）4.1多项式插值4.1.2插值多项式的存在唯一性因为

互不相同，所以式（4.2）不为零，根据解线性方程组的克莱姆(Cramer)法则，方程组的解

存在且唯一，从而

被唯一确定，这就证得了

次代数插值问题的解是存在且唯一的。（4.2）4.1多项式插值4.1.2插值多项式的存在唯一性基函数法：由线性空间的基出发，构造满足插值条件的多项式方法。用基函数法求插值多项式分两步：(1)定义

个线性无关的特殊代数多项式（插值基函数）

用

表示;(2)利用插值条件，确定插值基函数的线性组合表示的

次插值多项式：的系数 .（4.3）4.1多项式插值4.1.3