数值分析分解法
数值分析课程课件 直接三角分解方法
u22
u11
u2n
l n1 l n2
1
unn
即
a11 a12 a 21 a22
a1n
a2n
u11 l21u11
u12 l21u12 u22
u1n
l21u1n
u2n
a n1 a n2
ann
ln1u 11
由（5.3.1）- （5.3.4）求得L和U后，解方程组Ax=b 化为求解LUx=b，若记Ux=y，则有Ly=b。于是可分两部解 方程组LUx=b，只要逐次向前代入的方法即可求得y。第
二步求解Ux=y，只要逐次用向后回代的方法即可求得x。 设 x=（x1 ，x2， ···xn) T, y=（y1, y2, ···yn) T,
n
i1
lniuin
unn
第四章方程组的直接解法
由A的第1行和第1列可计算出U的第1行和L的第1列，即
u1 j a1 j , j 1, 2, , n,
（5.3.1）
lk1
ak1 u11
,k
2, 3,
, n.
（5.3.2）
如果U的第1至k-1列和L的第1至k-1列已经算出，则由
解 设 A=LU，即
l11 a11 1, l21 a21 2, l31 a31 0
u12
a12 l11
2, u13
a13 l11
1,
l22 a22 l21u12 3, l32 a32 l31u12 1
数值分析(06)矩阵分解法
(1) a11
( i k 1, , n) ( i k 1, , n;
(1) a12 ( 2) a 22 (1) a13 ( 2) a 23 ( 3) a 33
数值分析
LU分解的MATLAB程序 function A=lud(A) i k 1, , n %功能：对方阵A作三角分解A=LU,其中， (k ) (k ) % L为单位下三角阵， （1）m / akk U为上三角阵， ik aik %输入：方阵 A 。 (k 1) (k ) (k ) （ 2）aij aij mik akj ，j k 1, % 输出：紧凑存储 A=[L\U]. [n,n]=size(A); % 确定A的维数 for k=1:n-1 if A(k,k) ==0 break; end for i=k+1:n A(i,k) =A(i,k)/ A(k,k); A(i,k+1:n)= A(i,k+1:n)- A(i,k) *A(k,k+1:n); end end
1
(1) (1) a a 11 12 0 a ( 2) 22 0 a ( 2) n2 1 (1) a1, n ( 2) a2, n ( 3) a3,n A( 3) mi 2 ( 3) an ,n
(2)
... ... ...
其中分别是单位下上三角阵是对角定理矩阵分解定理数值分析数值分析非奇异知均非奇异而左边是单位下三角阵同时右边是对角阵故只能是单位阵而左边是下三角阵同时右边是单位上三角阵故只能是单位阵数值分析数值分析为对称正定阵则可唯一地分定理其中称正定分解三角阵
数值分析(07)矩阵的正交分解
H
(2) k
x
(2)
e
(k ) k 1
( 其中e1k ) (1, 0, , 0)T R n k 1 , 用前面介绍的方法 (2) 构造H k 。
U x y x i ei ( x1 , , xi i , , xn )T ,
有Hx y i ei
1 sign( x i ) 1
xi 0 xi 0
构造初等反射阵 UU T 1 T H I 2WW I 2 I UU T 2 U
解 : 3 sign( x3 ) x
2
4 0 4 1 3,因x3 2 0,
故取K 3 3 于是y 3e3 Ke3 (0, 0, 3, 0)T ,
U x y (2, 0, 5,1) , 3 ( 3 x3 ) 3(3 2) 15
数值分析
数值分析
function [H,y]=holder1(x) n=length(x); if x(1)<0 M=max(abs(x)); s=-s; if M==0, end; disp('M=0'); x(1)=s+x(1); return; p=s*x(1); else u=x; x=x/M; H=eye(n,n)-p\u*u'; end; y=zeros(n,1); s=norm(x); y(1)=-M*s;
1 T 1 2 2 其中 U U ( x1 ... ( xi i )2 xn ) 2 2 1 (2 xi i 2 i 2 ) i ( xi i ) 2
数值分析高斯顺序消去法、列主元消去法LU分解法
数值分析实验报告（1）学院：信息学院班级：计算机0903班姓名：***学号：********课题一A.问题提出给定下列几个不同类型的线性方程组，请用适当的方法求解线性方程组1、设线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------------------1368243810041202913726422123417911101610352431205362177586832337616244911315120130123122400105635680000121324⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125 x *= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 )T2、设对称正定阵系数阵线方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------------19243360021411035204111443343104221812334161206538114140231212200420424⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡87654321x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---4515229232060 x * = ( 1, -1, 0, 2, 1, -1, 0, 2 )T3、三对角形线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------4100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----5541412621357 x *= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )TB.（1）对上述三个方程组分别用Gauss 顺序消去法与Gauss 列主元消去法；平方根 与改进平方根法；追赶法求解（选择其一） （2）编写算法通用程序（3）在应用Gauss 消去时，尽可能利用相应程序输出系数矩阵的三角分解式C.(1)通过该课题的程序编制，掌握模块化结构程序设计方法 (2)掌握求解各类线性方程组的直接方法，了解各种方法的特点 (3)体会高斯消去法选主元的必要性 实验步骤：（高斯消去法，列主元，LU ）1顺序高斯消去法2.LU 分解法3.列主元高斯消去法（如下图）（1）高斯消去法运行结果如下（2）对方程的系数矩阵进行LU分解并求出方程组的解（3）列主元高斯消去法实验体会总结：利用gauss消去法解线性方程组的时候，如果没有经过选主元，可能会出现数值不稳定的现象，使得方程组的解偏离精确解。
数值分析实验报告
%消元过程
fori=k+1:n
m=A(i,k)/A(k,k);
forj=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<1e-10
flag='failure';return;
*x=(x0,x1….,xn),插值节点
*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值
*t求插值函数Pn(x)在t处的函数值
*返回值 插值函数Pn(x)在t处的函数值
*/
procedureNewton
forj=0to n
d1jyj;
endfor
forj=1to n
fori=j to n
[n,m]=size(A);nb=length(b)
%当方程组行与列的维数不相等时，停止计算，并输出出错信息
ifn~=m
error('The row and columns of matrix A must beepual!');
return;
end
%当方程组与右端项的维数不匹配时，停止计算，并输出错误信息
clear
fprintf('gauss-seidel迭代法')
x1_(1)=0;
x2_(1)=0;
x3_(1)=0;
fori=1:9
x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);
数值分析5LU分解法
数值分析5LU分解法LU分解法是一种常用的数值分析方法，用于解线性方程组。
本文将详细介绍LU分解法的原理、算法步骤、优缺点以及应用领域，以期能够全面地掌握这一方法。
一、LU分解法原理LU分解法是将一个方程组的系数矩阵A分解为两个矩阵L和U的乘积的形式，其中L是下三角矩阵，U是上三角矩阵，通过分解可以简化方程组的求解过程。
LU分解法的基本思想是将原始方程组Ax=b分解为Ly=b和Ux=y两个方程组，其中L和U是通过A分解得到的矩阵。
二、算法步骤1.首先，将系数矩阵A分解为两个矩阵L和U。
L是下三角矩阵，主对角线元素均为1，而U是上三角矩阵。
2.然后，将原始方程组Ax=b转化为Ly=b，求解y的值。
3.最后，将解y代入Ux=y，求解x的值，即可得到方程组的解。
三、算法优缺点1.优点：LU分解法将原始方程组的系数矩阵分解为两个形式简单的矩阵，简化了方程组的求解过程。
对于重复使用系数矩阵A的情况，只需要进行一次LU分解，然后根据新的b值求解新方程组，提高了计算效率。
2.缺点：LU分解法需要进行矩阵分解计算，计算量较大，因此对于规模较大的方程组计算效率较低。
此外，当系数矩阵A存在奇异性或病态时，LU分解法可能会失败。
四、应用领域LU分解法在科学计算领域有着广泛的应用，特别是在求解线性方程组方面。
例如，在工程领域中，常需要通过数值方法求解复杂的结构力学问题，此时可以使用LU分解法求解由有限元方法离散得到的大规模线性方程组。
另外，LU分解法还可以用于解非线性方程组、求逆矩阵、计算矩阵的行列式等。
总结：LU分解法是一种常用的数值分析方法，用于求解线性方程组。
通过将系数矩阵A分解为两个矩阵L和U的乘积形式，可以简化方程组的求解过程。
LU分解法的优点是提高了方程组的求解效率，适用于重复使用系数矩阵A的情况。
然而，LU分解法也存在一定的缺点，如计算量较大、对奇异性和病态问题的处理较为困难。
LU分解法在科学计算领域有着广泛的应用，可以用于求解工程问题中的大规模线性方程组，解非线性方程组，求逆矩阵等。
数值分析(计算方法)总结
第一章 绪论误差来源：模型误差、观测误差、截断误差（方法误差）、舍入误差ε(x )=|x −x ∗|是x ∗的绝对误差，e =x ∗−x 是x ∗的误差，ε(x )=|x −x ∗|≤ε，ε为x ∗的绝对误差限（或误差限） e r =ex =x ∗−x x为x ∗ 的相对误差，当|e r |较小时，令 e r =ex ∗=x ∗−x x ∗相对误差绝对值得上限称为相对误差限记为：εr 即：|e r |=|x ∗−x||x ∗|≤ε|x ∗|=εr绝对误差有量纲，而相对误差无量纲若近似值x ∗的绝对误差限为某一位上的半个单位，且该位直到x ∗的第一位非零数字共有n 位，则称近似值 x ∗有n 位有效数字，或说 x ∗精确到该位。
例：设x=π=3.1415926…那么x ∗=3,ε1(x )=0.1415926…≤0.5×100,则x ∗有效数字为1位，即个位上的3，或说 x ∗精确到个位。
科学计数法：记x ∗=±0.a 1a 2⋯a n ×10m (其中a 1≠0),若|x −x ∗|≤0.5×10m−n ，则x ∗有n 位有效数字，精确到10m−n 。
由有效数字求相对误差限：设近似值x ∗=±0.a 1a 2⋯a n ×10m （a 1≠0）有n 位有效数字，则其相对误差限为12a 1×101−n由相对误差限求有效数字：设近似值x ∗=±0.a 1a 2⋯a n ×10m （a 1≠0）的相对误差限为为12（a 1+1）×101−n 则它有n 位有效数字令x ∗、y ∗是x 、y 的近似值，且|x ∗−x|≤η(x )、|y ∗−y|≤η(y)1. x+y 近似值为x ∗+y ∗,且η(x +y )=η(x )+η（y ）和的误差（限）等于误差（限）的和2. x-y 近似值为x ∗−y ∗,且η(x +y )=η(x )+η（y ）3. xy 近似值为x ∗y ∗,η(xy )≈|x ∗|∗η(y )+|y ∗|∗η(x)4. η(xy )≈|x ∗|∗η(y )+|y ∗|∗η(x)|y ∗|21．避免两相近数相减2．避免用绝对值很小的数作除数 3．避免大数吃小数 4．尽量减少计算工作量 第二章 非线性方程求根1.逐步搜索法设f (a ) <0, f (b )> 0，有根区间为 (a , b )，从x 0=a 出发， 按某个预定步长(例如h =(b -a )/N )一步一步向右跨，每跨一步进行一次根的搜索，即判别f (x k )=f (a +kh )的符号，若f (x k )>0(而f (x k -1)<0),则有根区间缩小为[x k -1,x k ] (若f (x k )=0，x k 即为所求根), 然后从x k -1出发，把搜索步长再缩小，重复上面步骤，直到满足精度：|x k -x k -1|< 为止，此时取x *≈(x k +x k -1)/2作为近似根。
数值分析实验报告---高斯消去法 LU分解法
数值分析实验报告---高斯消去法 LU分解法实验一：高斯消去法一、实验目的1. 掌握高斯消去法的原理2. 用高斯消去法解线性方程组3. 分析误差二、实验原理高斯消去法(又称为高斯-约旦消去法)是一种利用矩阵消元的方法，将线性方程组化为改进的阶梯形式，从而解出线性方程组的解的方法。
具体而言，高斯消去法将线性方程组的系数矩阵化为一个上三角矩阵，再利用回带法求解线性方程组的解。
三、实验内容1.1、用高斯消去法解线性方程组在具体实验中，我们将使用高斯消去法来解决下述的线性方程组。
5x+2y+z=102x+6y+2z=14x-y+10z=25为了使用高斯消去法来解这个方程组，首先需要将系数矩阵A进行变换，消除A矩阵中第一列中的下角元素，如下所示:1, 2/5, 1/50, 28/5, 18/50, 0, 49/28接着使用回代法来计算该方程组的解。
回代法的过程是从下往上进行的，具体步骤如下：第三个方程的解：z=49/28；第二个方程的解： y=（14-2z-2x）/6；第一个方程的解： x=（10-2y-z）/5。
1.2、分析误差在使用高斯消去法求解线性方程组时，一般会出现截断误差，导致得到的解与真实解之间存在一些误差。
截断误差的大小和矩阵的维数有关。
为了估计截断误差，我们使用矩阵B来生成误差，在具体实验中，我们将使用下面的矩阵：我们来计算该矩阵的行列式，如果方程组有唯一解，则行列性不为0。
本例中，行列式的值是 -1，因此方程组有唯一解。
然后我们计算真实解和高斯消去法得到的解之间的误差，具体公式如下所示：误差 = 真实解的范数 - 高斯消去法得到的解的范数其中，范数的定义如下：||x||1=max{|xi|}; ||x||2=sqrt{(|x1|^2 + |x2|^2 + ... + |xn|^2)}四、实验步骤1、将高斯消去法的每一个步骤翻译成代码，并保存为一个独立的函数。
2、将代码上传至 Python 交互式环境，并使用高斯消去法来解线性方程组。
(完整版)数值分析重点公式
第一章 非线性方程和方程组的数值解法 1）二分法的基本原理，误差：~12k b ax α+--<2）迭代法收敛阶：1lim0i pi ic εε+→∞=≠，若1p =则要求01c <<3）单点迭代收敛定理：定理一：若当[],x a b ∈时，[](),x a b ϕ∈且'()1x l ϕ≤<，[],x a b ∀∈，则迭代格式收敛于唯一的根；定理二：设()x ϕ满足：①[],x a b ∈时，[](),x a b ϕ∈， ②[]121212,,, ()(),01x x a b x x l x x l ϕϕ∀∈-≤-<<有 则对任意初值[]0,x a b ∈迭代收敛，且：110111i i iii x x x ll x x x lαα+-≤---≤--定理三：设()x ϕ在α的邻域内具有连续的一阶导数，且'()1ϕα<，则迭代格式具有局部收敛性；定理四：假设()x ϕ在根α的邻域内充分可导，则迭代格式1()i i x x ϕ+=是P 阶收敛的 ()()()0,1,,1,()0j P j P ϕαϕα==-≠L （Taylor 展开证明）4）Newton 迭代法：1'()()i i i i f x x x f x +=-，平方收敛 5）Newton 迭代法收敛定理：设()f x 在有根区间[],a b 上有二阶导数，且满足： ①：()()0f a f b <； ②：[]'()0,,f x x a b ≠∈；③：[]'',,f x a b ∈不变号④：初值[]0,x a b ∈使得''()()0f x f x <；则Newton 迭代法收敛于根α。
6）多点迭代法：1111111()()()()()()()()()i i i i i i i i i i i i i i i f x f x f x x x x x f x f x f x f x f x f x x x -+-----=-=+----收敛阶：P =7）Newton 迭代法求重根（收敛仍为线性收敛），对Newton 法进行修改 ①：已知根的重数r ，1'()()i i i i f x x x rf x +=-（平方收敛） ②：未知根的重数：1''()(),()()()i i i i u x f x x x u x u x f x +=-=，α为()f x 的重根，则α为()u x 的单根。
Doolittle分解
矩阵数值分析实验报告专业信息与计算科学班级学号姓名指导教师Doolittle 分解法一、实验目的在Gauss 消元法中,对于n 阶方程组,应用消去发经过n-1步消元之后,得到一个与Ax=b 等价的代数线性方程组)1()1(--=n n b x A ,而且)1(-n A 为一个上三角矩阵.所以我们想是否能把矩阵A 分解成一个下三角阵与一个上三角阵的乘积 A=LR,其中L 为下三角阵,R 为上三角阵.就变成了两个三角形方程组⎩⎨⎧==yRx b Ly , 的求解问题。
二、算法思想Setp1：利用for 循环求出r[k][j]=a[k][j]-∑-=1k 1p ]kp [r ]kp [l ，l[ik]=（a[ik]-∑-=1k 1p ]ip [r ]ip [l ）/r[k][k]。
Step2：y[i]=b[i]-∑-=1i 1k ]k [y ]ik [l ,得出x[i]=(y[i]-∑+=n1i k ]k [x ]ik [r )/r[i][i].三、程序代码#include <stdio.h>#include <stdlib.h>#define N 10 //矩阵大小范围float getmx(float a[N][N], float x[N], int i, int n) {float mx = 0;int r;for(r=i+1; r<n; r++){mx += a[i][r] * x[r];}return mx;}float getmy(float a[N][N], float y[N], int i, int n) {float my = 0;int r;for(r=0; r<n; r++){if(i != r) my += a[i][r] * y[r];}return my;}float getx(float a[N][N], float b[N], float x[N], int i, int n) {float result;if(i==n-1) //计算最后一个x的值result = (float)(b[i]/a[n-1][n-1]);else //计算其他x值（对于公式中的求和部分，需要调用getmx()函数） result = (float)((b[i]-getmx(a,x,i,n))/a[i][i]);return result;}float gety(float a[N][N], float b[N], float y[N], int i, int n) {float result;if(i==0) //计算第一个y的值result = float(b[i]/a[i][i]);else //计算其他y值（对于公式中的求和部分，需要调用getmy()函数） result = float((b[i]-getmy(a,y,i,n))/a[i][i]);return result;}void main(){float l[N][N]={0}; //定义L矩阵float u[N][N]={0}; //定义U矩阵float y[N]={0}; //定义数组Yfloat x[N]={0}; //定义数组Xfloat a[N][N]={{0},{0},{0}}; //定义系数矩阵float b[N]={0}; //定义右端项float sum=0;int i,j,k;int n=0;int flag=1;//用户手工输入矩阵while(flag){printf("请输入系数矩阵的大小：");scanf("%d", &n);if(n>N){printf("矩阵过大！\n");continue;}flag=0;}printf("请输入系数矩阵值：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("a[%d][%d]: ", i, j);scanf("%f", &a[i][j]);}}printf("请输入右端项数组：\n");for(i=0; i<n; i++){printf("b[%d]: ", i);scanf("%f", &b[i]);}/*显示原始矩阵*/printf("原始矩阵：\n");for(i=0; i<n; i++){for(j=0; j<n; j++)printf("%0.3f ",a[i][j]);printf("\n");}printf("\n\n");/*初始化矩阵l*/for(i=0; i<n; i++){for(j=0; j<n; j++){if(i==j) l[i][j] = 1;}}/*开始LU分解*//*第一步：对矩阵U的首行进行计算*/for(i=0; i<n; i++){u[0][i] = (float)(a[0][i]/l[0][0]); }/*第二步：逐步进行LU分解*/for(i=0; i<n; i++){/*对“L列”进行计算*/for(j=i+1; j<n; j++){for(k=0,sum=0; k<n; k++){if(k != i) sum += l[j][k]*u[k][i];}l[j][i] = (float)((a[j][i]-sum)/u[i][i]); }/*对“U行”进行计算*/for(j=i+1; j<n; j++){for(k=0,sum=0; k<n; k++){if(k != i+1) sum += l[i+1][k]*u[k][j]; }u[i+1][j] = (float)((a[i+1][j]-sum));}}/*输出矩阵l*/printf("矩阵L：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%0.3f ", l[i][j]);}printf("\n");}/*输出矩阵u*/printf("\n矩阵U：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%0.3f ", u[i][j]);}printf("\n");}/*回代方式计算数组Y*/for(i=0; i<n; i++){y[i] = gety(l,b,y,i,n);}/*显示数组Y*/printf("\n\n数组Y：\n");for(i=0; i<n; i++){printf("y%d = %0.3f\n", i+1,y[i]); }/*回代方式计算数组X*/for(i=n-1; i>=0; i--){x[i] = getx(u,y,x,i,n);}/*显示数组X*/printf("\n\n数组X：\n");for(i=0; i<n; i++){printf("x%d = %0.3f\n", i+1,x[i]); }}四、运行结果五、参考文献[1]刑志栋，矩阵数值分析，陕西：陕西科学技术出版社， 2005。
研究生数值分析(8)分解
b1 bk
k 1
lkj
yj
j 1
(k 2,3,
（5）
, n)
xn
xk
yn / unn (yk
n
ukj xj ) / ukk
j k 1
（6）
(k n 1,n 2, ,1)
杜利特尔矩阵分解 求解线性方程组的过程为: 10 实现A=LU分解，即
(a)按计算公式(1),(2)依次计算U的第1行元素
(i k, k 1, , n) （3）
k 1
lik (aik liju jk ) / ukk j 1
(i k 1, , n; k n) （4）
在我们利用杜利特尔矩阵分解解线性方程 组AX=b时,只要实现矩阵分解A=LU,依次解三角 形方程组LY=b与UX=Y即可。
计算公式：
y1
yk
u1i (i 1, 2, , n) 与L的第1列元素 li1 (i 2,3, , n)
(b) 对k+2,3,…,n 按计算公式(3),(4)依次
计算U的第k行元素 uki (i k, k 1, , n) 与L的第
k列元素 lik (i k 1, , n; k n)
20 求解三角形方程组LY=b，即按计算公
若U为单位上三角阵(对角元都是1的上三角阵),
L为下三角阵,则称为克劳特(Crout)分解。
下面分析实现矩阵杜利特尔(Doolittle)分解和克 劳特(Crout)分解的条件,讨论这些分解的唯一性。
定理2 (矩阵三角分解基本定理)
设 A Rnn 。若A的顺序主子式
det( Ak ) 0
(k 1, 2, , n)
li,k 1
lik
ln,k 1
数值分析LU分解法
数值分析LU分解法数值分析是计算机科学中的一个重要分支，它研究数值计算的原理、方法和算法。
数值分析的一个常见问题是解线性方程组，而LU分解法是一种常用的解线性方程组的数值方法。
本文将详细介绍LU分解法的原理、应用和优势。
LU分解法是一种矩阵分解的方法，它将一个矩阵分解为两个特定的矩阵：L下三角矩阵和U上三角矩阵。
具体而言，给定一个n×n的矩阵A，LU分解法可以将其表示为LU=A，其中L是一个n×n的下三角矩阵，U是一个n×n的上三角矩阵。
通过将矩阵A分解为L和U的乘积形式，LU分解法可以更高效地解线性方程组。
LU分解的过程大致如下：1.初始化L和U为n×n的零矩阵。
2.对于矩阵A的第一列，将其元素依次作为L的第一列。
3.对于矩阵A的第一行，将其元素依次作为U的第一行。
4.对于矩阵A的每一列j（除了第一列），计算L矩阵的第j列和U矩阵的第j行。
具体方法是通过消元法将L矩阵的第j列元素计算为A矩阵的第j列元素除以U矩阵的第j行元素。
LU分解法的主要优势在于解线性方程组的效率和稳定性。
与高斯消元法相比，LU分解法具有更好的稳定性，特别是在矩阵的行主元不为零时。
当需要多次求解不同右端向量的线性方程组时，通过LU分解可只需一次分解，多次使用L和U来求解不同的线性方程组。
这样可以大大提高计算效率。
此外，LU分解法还有其他一些应用。
例如，通过LU分解法可以计算矩阵的逆和行列式。
由于LU分解法将矩阵分解为两个三角矩阵，求解矩阵的逆只需分别求解L和U的逆并将其乘积即可。
同样地，通过LU分解法也可以计算矩阵的行列式，行列式的值等于U矩阵对角线上的元素的乘积。
总结来说，LU分解法是一种在数值分析中常用的解线性方程组的方法。
它通过将矩阵分解为下三角矩阵和上三角矩阵的乘积形式，实现了高效和稳定的线性方程组求解。
此外，LU分解法还可以用于计算矩阵的逆和行列式。
因此，LU分解法在数值计算和科学计算中有着广泛的应用。
数值分析-第二章小结
第二章 线性方程组的数值解法-------学习小结姓名 班级 学号 一、本章学习体会通过本章的学习，我了解了线性方程组的不同解法，切实体会到了不同的计算方法对计算结果的影响。
求解线性方程组的方法可分为两大类：直接方法和迭代方法。
直接方法在解一般的线性方程组的时候比较简便，使用此方法经过有限次运算就可得到方程组的解。
然而迭代法是要构造一个无限的向量序列，其极限是方程组的解向量，它适用于求解大型稀疏线性方程组。
总的来说，直接方法和迭代法各有优点与不足，在解线性方程组的时候，我们要根据具体的线性方程组的特点来选择合适的解法，这样我们才能快速准确的得到方程组的解。
因此，我们要熟悉书中介绍的各类线性方程组的解法，同时要善于思考、总结，在使用各种方法求解的同时尽量提出自己独特的见解，通过不断练习计算，使自己的能力得到提高。
二、本章知识梳理线性方程组的求解方法分为直接法和迭代法两种，Gramer （克莱姆）法是直接法的一种，但由于其计算量比较大，在世界工作中其效率比较低、经济效益差，所以此方法我们很少使用，本章主要介绍其他的计算方法。
2.1 Gauss 消去法Gauss （高斯）消去法由消元和回代两个过程组成。
消元过程就是对方程组的增广矩阵做有限次的初等行变换，使它的系数矩阵部分变换为上三角阵。
所用的初等行变换主要有两种：第一种，交换两行的位置；第二种，用一个数乘某一行加到另一行上。
回代过程就是先由方程组的最后一个方程解出n x ，然后通过逐步回代，依次求出1n x -，2n x -，…，1x 。
这种Gauss 消去法可分为Gauss 消去法和列主元素Gauss 消去法两种。
2.1.1 顺序Gauss 消去法在Gauss 消去法的消元过程中对方程组的增广矩阵只做前述的第二种初等行变换就形成了顺序Gauss 消去法，其算法如下：记(1)ij ij a a = （i ，j=1，2，…，n ）i i 1、 消元过程对于k=1,2，…，n-1执行 （1）如果()0k kka =，则算法失效，停止计算；否则转（2）。
数值分析底部剪力法与振型分解反应谱法对比分析
底部剪力法与振型分解反应谱法对比分析摘要：建筑结构抗震设计是建筑结构设计中必不可少，也是非常重要的一部分。
结构抗震在建筑结构的总成本中占有相当大的比例。
建筑抗震设计规范中有关于结构抗震计算的方法以及适用范围，水平地震力的计算方法主要是底部剪力法和振型分解反应谱法，底部剪力法适用于质量和刚度沿高度分布比较均匀的结构，而振型分解反应谱法能反应结构的真实情况，对一般结构都适用。
本文通过对五层、八层、十层，质量和刚度分布均匀和不均匀框架结构的各层剪力计算，来比较两种方法的计算结果，验证底部剪力法的适用范围以及有效性。
本文对结构特征周期的计算是用广义Jacobi方法，通过Fortran语言编程实现的。
关键词：底部剪力法；振型分解反应谱法；Jacobi方法；Fortran语言Comparative Analysis between Equivalent Base Shear Method and ModalAnalysis MethodAbstract: Seismic design plays an essential and important part in the structure design. It also makes up a significant proportion of the total cost. About the horizontal seismic force, the code has detailed specification of the calculation principle and applicable scope. The calculation method for horizontal seismic force mainly has the equivalent base shear method and modal analysis method. The equivalent base shear method is suitable for mass and stiffness along the height of structure with uniform distribution, and the modal analysis method reflects the true action of the structure and has a wide usage. By calculating the shear of five-story, eight-story and ten-story framework with mass uniform or non-uniform distribution, this paper verified the scope and the effectiveness of the equivalent base shear method. The eigenperiod of the structure is calculated by generalized Jacobi method though the Fortran language programming.Key words: Equivalent Base Shear Method; Modal Analysis Method; Jacobi Method; Fortran Language引言实际的建筑结构其质量一般均是连续分布的，因此，严格的说，其动力自由度均是无限的。
数值分析公式大全
的位置未知，但有截断误差限：
(a,b)
3， 均差（差商） 一阶均差；f[x0，xk]=
,Mn+1=
二阶均差：f[x0,，x1，xk]= ，
，
高阶均差：f[x0,，x1，…，xk]= ， ， ，
，， ，
性质：1，k 阶均差可表示为函数值 f（x0），f（x1），…，f（xn）的线性组合 2，对称性，与节点次序无关
Ax=b 将 A 按行化简为三角矩阵（等同于做多次消元过程）最后解简单方程组 A（n）x=b（n） 2， 高斯主元素消去法 列主元素消去法：若出现 akk（k）=0 B= 在 A 的第一列中选择绝对值最大元素做为主元素，如丨 ai1，1 丨=max1≤i≤n 丨 ai1 丨然后 交换 B 的第一行与第 i1 行， →
≤
5， 差分
等距离节点 xk=x0+kh，k=0，1，…，n；fk=f（xk）
xk 处的一阶向前差分：Δ fk=fk-+1-fk，xk 处的二阶向前差分：Δ 2fk=Δ fk-+1-Δ fk；
xk 处的
n
阶差分：Δ
nfk=Δ
n-1fk-+1-Δ
f n-1 k
【差分与差商的关系】f[xk,xk+1]=
余项 Rn= 4，高斯-勒让德求积公式
程组中可求
，其中高斯点为 Pn+1（x）=0 的解，将 代入高斯公式所得的方
Rn= 5，高斯-切比雪夫求积公式
，其中高斯点为 Tn+1（x）=0 的解，
k=0，1，…，n。Ak=
也可写为
，
，k=1，2，…，n
第五章解线性方程 组的直接方法：
去除矩阵论部分的 基本知识点，剩余内容有 ; 1， 高斯消去法
清华大学研究生高等数值分析计算实验奇异值分解SVD以及图像压缩matlab源程序代码
第1部分方法介绍奇异值分解（SVD ）定理：设m n A R ⨯∈，则存在正交矩阵m m V R ⨯∈和n n U R ⨯∈，使得TO A V U O O ∑⎡⎤=⎢⎥⎣⎦其中12(,,,)r diag σσσ∑= ，而且120r σσσ≥≥≥> ，(1,2,,)i i r σ= 称为A 的奇异值，V 的第i 列称为A 的左奇异向量，U 的第i 列称为A 的右奇异向量。
注：不失一般性，可以假设m n ≥，（对于m n <的情况，可以先对A 转置，然后进行SVD 分解，最后对所得的SVD 分解式进行转置，就可以得到原来的SVD 分解式）方法1：传统的SVD 算法主要思想：设()m n A R m n ⨯∈≥，先将A 二对角化，即构造正交矩阵1U 和1V 使得110T B n U AV m n⎡⎤=⎢⎥-⎣⎦其中1200n n B δγγδ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦然后，对三角矩阵T T B B =进行带Wilkinson 位移的对称QR 迭代得到:T B P BQ =。
当某个0i γ=时，B 具有形状12B O B O B ⎡⎤=⎢⎥⎣⎦，此时可以将B 的奇异值问题分解为两个低阶二对角阵的奇异值分解问题；而当某个0i δ=时，可以适当选取'Given s 变换，使得第i 行元素全为零的二对角阵，因此，此时也可以将B 约化为两个低阶二对角阵的奇异值分解问题。
在实际计算时，当i B δε∞≤或者()1j j j γεδδ-≤+（这里ε是一个略大于机器精度的正数）时，就将i δ或者i γ视作零，就可以将B 分解为两个低阶二对角阵的奇异值分解问题。
主要步骤：（1）输入()m n A R m n ⨯∈≥及允许误差ε（2）计算Householder 变换1,,,n P P ⋅⋅⋅，12,,n H H -⋅⋅⋅，使得112()()0Tn n B nP P A H H m n -⎡⎤⋅⋅⋅⋅⋅⋅=⎢⎥-⎣⎦其中1200n n B δγγδ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ 12:n U PP P =⋅⋅⋅，122:.n V H H H -=⋅⋅⋅ （3）收敛性检验：（i ）将所有满足()1j j j γεδδ-≤+的j γ置零；（ii ）如果0,2,,j j n γ== ，则输出有关信息结束；否则，1:0γ=，确定正整数p q <，使得10p q n γγγ+==⋅⋅⋅==，0j γ≠，p j q <≤；（iii ）如果存在i 满足1p i q ≤≤-使得i B δε∞≤，则:0i δ=，1:i x γ+=，1:i y δ+=，1:0i γ+=，:1l =，转步（iv ）；否则转步（4） （iv ）确定cos ,sin c s θθ==和σ使0c s x s c y σ⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦这也相应于0Tc s y s c x σ⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦所以可以直接调用'Given s 变换算法得到:i l δσ+=，:(,,)T U UG i i l θ=+这相当于(1:;,)(1:;,)(1:;,)Tc s c s U n i i l U n i i l U n i i l s c s c -⎡⎤⎡⎤+=+=+⎢⎥⎢⎥-⎣⎦⎣⎦（v ）如果l q i <-，则1:i l x s γ++=，11:i l i l c γγ++++=，1:i l y δ++=，:1l l =+转步（iv ），否则转步（i ）（4）构造正交阵P 和Q ，使12=T P B Q B 仍为上双对角阵，其中112100pp p p q q B δγδγγδ+++⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦， 得121:=T B B P B Q =，:(,,)p n p q U Udiag I P I --=，:(,,)p n p q V Vdiag I Q I --=然后转步（3）。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
§3 LU 分解法——Gauss 消去法的变形知识预备：1矩阵的初等行变换、初等矩阵及其逆、乘积2矩阵的乘法3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆一、Gauss 消去法的矩阵解释Gauss 消去法实质上是将矩阵A 分解为两个三角矩阵相乘。
我们知道，矩阵的初等行变换实质就是左乘初等矩阵。
第一轮消元：相当于对A（1）左乘矩阵L 1，即)2()1(1A A L =其中)1(11)1(11)2()2(2)2(2)2(22)1(1)1(12)1(11)2(131211,,1001011a a l a a a a a a a Al l l L i i nn n n n n =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=第二轮消元：对应于)3()2(2A A L =一般地1,,2,1)1()(-==+n k A A L k k k (1)其中nk k i a a l l l L k kkk ikik nk kk k ,,2,1,,10111)()(1 ++==⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=+整个消元过程为U A A L L L L n n n 记)(1221=-- ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=nn nnu u u u u u 22211211………（2） 从而U L U L L L L U L L L L A n n n n ⋅===---------1112121111221)(其中L 是单位下三角矩阵，即,1,,1,,3,2,,1111)()(21323121⎪⎪⎭⎫⎝⎛-===⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=n j n i l l l l l l L j jj j ija a ij n n …（3） 【注】消元过程等价于A 分解成LU 的过程回代过程是解上三角方程组的过程。
二、矩阵的三角分解1、若将A 分解成L?U ，即A=L?U ，其中L 为单位下三角矩阵,U 为非奇异上三角矩阵,则称之为对A 的Doolittle 分解。
当A 的顺序主子式都不为零时，消元运算可进行，从而A 存在唯一的Doolittle 分解。
证明：若有两种分解，A=L 1U 1，A=L 2U 2，则必有L 1=L 2，U 1=U 2。
因为L 1U 1=L 2U 2，而且L 1，L 2都是单位下三角矩阵，U 1,U 2都是可逆上三角矩阵,所以有112112--=U U L L因此 )(112112单位矩阵I U U L L ==--即L 1=L 2，U 1=U 2、2、若L 是非奇异下三角矩阵，U 是单位上三角矩阵时，A 存在唯一的三角分解，A=LU ，称其为A 的Crout 分解（对应于用列变换实施消元）三、直接分解（LU 分解）算法LU 分解算法公式——按矩阵乘法⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=nn nnn n u u u u u u l l l l l A 22211211213231211111 第一步：利用A 中第一行、第一列元素确定U 的第一行、L 的第一列元素。
由),,2,1()0,,0,,,()0,,0,0,1(1211n j u u u u a jT ij j j j ==⋅=),,3,2()0,,0,()0,,1,,,(111111211n i u l u l l l a i T ii i i i =⋅=⋅=-得 u 1j =a 1j n j ,,2,1( =) l i1=a i1/u 11n i ,,3,2( =)第r 步:利用A 中第r 行、第r 列剩下的元素确定U 的第r 行、L 的第r 列元素（r=2，3，…，n ）.由),,1,()0,,0,,,()0,,0,1,,,(1121121n r r j u u l u u u l l l a rjkj r k rk Tjj j j rr r r rj +=+=⋅=∑-=-得U 的第r 行元素为∑-=+=-=11,,1,,r k kj rk rj rj n r r j u l a u由),,2,1()0,,0,,,()0,,0,1,,,(1121121n r r i u l u l u u u l l l a rrir kr r k ik Trr r r ii i i ir ++=+=⋅=∑-=-得),,1,1,,3,2(/)(11n r i n r u u l a l rrkr r k ik ir ir +=-=-=∑-=…………(4)直接分解的紧凑格式：方程组的三角分解算法（LU 分解）对于方程组Ax=b ，设A=LU (Doolittle 分解)。
由于 ⎩⎨⎧==⇔=yUx bLy b Ax1、求解Ly=b ：∑-==-==1111),,3,2(,,i k k ik i i n i y l b y b y (5)2、求解Ux=y ：)1,2,,2,1(,/)(,/1--=-==∑+=n n i u x uy x u y x ii ni k k iki i nn n n (6)LU 分解算法步1，输入A ，b ；步2，对j=1,2,…,n 求,:111j j j a u u = 对i=2,3,…,n 求;/:11111u a l l i i i = 步3,对r=2,3,…,n 做（）-: （）∑-=+=-=11),,,1,(,r k kj rkrj rj n r r j u la u（）);;,,1(/)(11n r n r i u u la l rrr k kr ikir ir ≠+=-=∑-=步4，∑-=-===1111;:,,,3,2,i k k iki i i y lb y y n i b y 求对步5，ii ni k k iki i i n n u x uy x x n i y x /)(:1,,1,1∑+=-=-==求对步6，输出);,,2,1(n i x i =结束。
⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-713542774322322x x x 解：对系数矩阵A 进行LU 分解6)(,2,2/)(1,3,1,2,3,2,22332133133332222123132322322121223121131211=--===-===-=-=====u l u l a u u u u l a l u u u l a u l l u u u j j j 所以因此⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=613322121121A 先解627,521,3,213121=-+-=-=-===y y y y y y b Ly 则。
再解22/)323(,23/)5(,1,321323=--=-=--===x x x x x x y Ux 解出程序：LU_factorization%Not Select Column LU_factorization clear alln=3;a=[2 2 3;4 7 7;-2 4 5];b=[3;1;-7]; %n=3;a=[1 4 7;2 5 8;3 6 11];b=[1;1;1]; %LU_factorazation for i=2:na(i,1)=a(i,1)/a(1,1); end afor r=2:n for j=r:n s=0.;for k=1:r-1s=s+a(r,k)*a(k,j); enda(r,j)=a(r,j)-s; endfor i=r+1:n s=0.;for k=1:r-1s=s+a(i,k)*a(k,r); enda(i,r)=(a(i,r)-s)/a(r,r); end a end%Extract Lower/Upper Triangular Part l=tril(a); for i=1:n l(i,i)=1; endu=triu(a); l u%Linear Lower Triangular Equation Solution y=l\b%Linear Upper Triangular Equation Solution x=u\y四、列主元LU 分解当用LU 分解法解方程组时，从第r(r=1,2,…,n)步分解计算公式可知 ∑-=+=-=11),,1,(r k kjrkrj rj n r r j u la u),,1(/)(11n r i u u l a l rrkr r k ik ir ir +=-=∑-=当rr u 很小时，可能引起舍入误差的累积、扩大。
因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列主元消去法在理论上是等价的），它通过交换A 的行实现三角分解PA=LU ，其中P 为置换阵。
设第r-1步分解计算己完成，则有⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡→---r n n r r r r n l l l u u u lu u A,1,,12221111第r 步计算时为了避免用绝对值很小的数作除数，引进中间量：∑-==-=11),,(,r k krik ir i n r i u l a S则有：),,1(,n r i S S l S u r i ir r rr +===（1） 选主元：确定i ni r ir r S S i ≤≤=m ax ,使（2） 交换两行：,时当r i r ≠交换A 的第r 行与第r i 行元素（相当于先交换原始矩阵A 第r 行与第r i 行元素后，再进行分解计算得到的结果，且1≤ir l ）（3） 进行分解计算 附：列主元LU 分解a=[2 2 3;4 7 7;-2 4 5];b=[3;1;-7]; [l,u]=lu(a) y=l\b x=u\y%x=inv(u)*inv(v)*b[l,u,p]=lu(a) y=l\(p*b) x=u\y。