LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序
一、实验目的及题目
1.1 实验目的：
（1）学会用高斯列主元消去法，LU 分解法，Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。
（2）学会用Matlab 编写各种方法求解线性方程组的程序。 1.2 实验题目：
1. 用列主元消去法解方程组：
1241234
123412343421233234
x x x x x x x x x x x x x x x ++=??+-+=??--+=-??-++-=?
2. 用LU 分解法解方程组,Ax b =其中
4824012242412120620266216A --?? ?- ?= ? ?-??，4422b ?? ? ?= ?- ?-??
3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组：
123234
1231234102118311210631125
x x x x x x x x x x x x x -+=-??-+=-??-+=??-+-+=?
二、实验原理、程序框图、程序代码等
2.1实验原理
2.1.1高斯列主元消去法的原理
Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：
1111221122222n n n n nn n n b x b x b x g b x b x g b x g +++=??++=?
?
?
?=?
这个过程就是消元，然后再回代就好了。具体过程如下：
对于1,2,,1k n =- ，若()
0,k kk a ≠依次计算
()()
(1)()()(1)()()/,,1,,k k ik ik kk k k k ij ij ik kj
k k k i i ik k m a a a a m a b b m b i j k n
++==-=-=+
然后将其回代得到：
()()
()()()1/()/,1,2,,1n n n n nn n k k k k
k kj j kk j k x b a x b a x a k n n =+?=??=-=--??
∑ 以上是高斯消去。
但是高斯消去法在消元的过程中有可能会出现()
0k kk a =的情况，这时消元就无法进行了，即使主元数()
0,k kk a ≠但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差
的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。 2.1.2直接三角分解法（LU 分解）的原理
先将矩阵A 直接分解为A LU =则求解方程组的问题就等价于求解两个三角形方程组。 直接利用矩阵乘法，得到矩阵的三角分解计算公式为：
1111111
1
1
1,1,2,,/,2,,,,,1,,,2,3,()/,1,2,,i i i i k kj kj km mj m k ik ik im mk kk
m u a i n
l a u i n
u a l u j k k n k n l a l u u i k k n k n
-=-===??
==??=-=+??=??=-=++≠??
∑∑ 且 由上面的式子得到矩阵A 的LU 分解后，求解Ux=y 的计算公式为
11
111,2,3,/()/,1,2,,1
i i i ij j j n n nn n i i ij j ii j i y b y b l y i n
x y u x y u x u i n n -==+=???
=-=??
=???
=-=--??
∑∑
以上为LU 分解法。
2.1.3Jacobi 迭代法和Gauss-Seidel 迭代法的原理 （1）Jcaobi 迭代
设线性方程组
b Ax = (1)
的系数矩阵A 可逆且主对角元素nn a ,...,a ,a 2211均不为零,令 ()nn a ,...,a ,a diag D 2211=
并将A 分解成
()D D A A +-= (2) 从而(1)可写成
()b x A D Dx +-= 令
11f x B x +=
其中b D f ,A D I B 1
111
--=-=. (3) 以1B 为迭代矩阵的迭代法(公式)
()()111f x B x k k +=+ (4)
称为雅可比(Jacobi)迭代法,其分量形式为
?
??
[]
,...,,k ,
n ,...,i x a b
a x
n
i
j j )
k (j j i i
ii
)k (i
21021111==∑-=≠=+ (5)
其中
()()()()
()T
n x ,...x ,x x 002010=为初始向量. （2）Gauss-Seidel 迭代
由雅可比迭代公式可知,在迭代的每一步计算过程中是用()
k x 的全部分量来计算()
1+k x
的
所有分量,显然在计算第i 个分量()
1+k i x 时,已经计算出的最新分量()
()
11
11
+-+k i k x ,...,x 没有被利
用。
把矩阵A 分解成
U L D A --= (6)
其中()nn
a ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下三角和上三角部分,
于是,方程组(1)便可以写成
()b Ux x L D +=- 即
22f x B x +=
其中
()()b L D f ,
U L D B 1
21
2---=-= (7)
以2B 为迭代矩阵构成的迭代法(公式)
()()221f x B x k k +=+ (8)
称为高斯—塞德尔迭代法,用分量表示的形式为
?
??[]
,...
,,k ,n ,,i x a x a b a x
i j n i j )
k (j ij )
k (j ij i ii
)k (i
21021111111==∑∑--=-=+=++
2.2程序代码
2.2.1高斯列主元的代码
function Gauss(A,b) %A 为系数矩阵，b 为右端项矩阵 [m,n]=size(A); n=length(b); for k=1:n-1
[pt,p]=max(abs(A(k:n,k))); %找出列中绝对值最大的数 p=p+k-1; if p>k
t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; %交换行使之变到主元位置上 t=b(k);b(k)=b(p);b(p)=t; end
m=A(k+1:n,k)/A(k,k); %开始消元 A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag~=0
Ab=[A,b];
end
end
x=zeros(n,1); %开始回代
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end
for k=1:n
fprintf('x[%d]=%f\n',k,x(k));
end
2.2.2 LU分解法的程序
function LU(A,b) %A为系数矩阵，b为右端项矩阵
[m,n]=size(A); %初始化矩阵A，b，L和U
n=length(b);
L=eye(n,n);
U=zeros(n,n);
U(1,1:n)=A(1,1:n); %开始进行LU分解
L(2:n,1)=A(2:n,1)/U(1,1);
for k=2:n
U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k); end
L %输出L矩阵
U %输出U矩阵
y=zeros(n,1); %开始解方程组Ux=y
y(1)=b(1);
for k=2:n
y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
end
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end
for k=1:n
fprintf('x[%d]=%f\n',k,x(k));
end
2.2.3 Jacobi迭代法的程序
function Jacobi(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度
[m,n]=size(A);
D=diag(diag(A)); %求矩阵D
L=tril(A)-D; %求矩阵L
U=triu(A)-D; %求矩阵U
temp=1;
x=zeros(m,1);
k=0;
while abs(max(x)-temp)>eps
temp=max(abs(x));
k=k+1; %记录循环次数
x=-inv(D)*(L+U)*x+inv(D)*b; %雅克比迭代公式
end
for k=1:n
fprintf('x[%d]=%f\n',k,x(k));
end
2.2.4Gauss-Seidel迭代程序
function Gauss_Seidel(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);
D=diag(diag(A)); %求矩阵D
L=D-tril(A); %求矩阵L
U=D-triu(A); %求矩阵U
temp=1;
x=zeros(m,1);
k=0;
while abs(max(x)-temp)>eps
temp=max(abs(x));
k=k+1; %记录循环次数
x=inv(D-L)*U*x+inv(D-L)*b; %Gauss_Seidel的迭代公式
end
for k=1:n
fprintf('x[%d]=%f\n',k,x(k));
end
三、实验过程中需要记录的实验数据表格
3.1第一题（高斯列主元消去）的数据
>> A=[1 1 0 3;2 1 -1 1; 3 -1 -1 3;-1 2 3 -1];
>> b=[4;1;-3;4];
>> Gauss(A,b)
x[1]=-1.333333
x[2]=2.333333
x[3]=-0.333333
x[4]=1.000000
3.2第二题（LU分解法）数据
>> A=[48 -24 0 -12;-24 24 12 12;0 6 20 2;-6 6 2 16];
>> b=[4; 4;-2;-2];
>> LU(A,b)
L =
1.0000 0 0 0
-0.5000 1.0000 0 0
0 0.5000 1.0000 0
-0.1250 0.2500 -0.0714 1.0000
U =
48.0000 -24.0000 0 -12.0000
0 12.0000 12.0000 6.0000
0 0 14.0000 -1.0000
0 0 0 12.9286
x[1]=0.521179
x[2]=1.005525
x[3]=-0.375691
x[4]=-0.259669
3.3第三题Jacobi迭代法的数据
>> A=[10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11];
b=[-11;-11;6;25];
Jacobi(A,b,0.00005)
x[1]=-1.467396
x[2]=-2.358678
x[3]=0.657604
x[4]=2.842397
3.4第三题用Gauss_Seidel迭代的数据
>> A=[10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11];
>> b=[-11;-11;6;25];
>> Gauss_Seidel(A,b,0.00005)
x[1]=-1.467357
x[2]=-2.358740
x[3]=0.657597
x[4]=2.842405
四、实验中存在的问题及解决方案
4.1存在的问题
（1）第一题中在matlab中输入>> Gauss(A,b)（数据省略）得到m =4 n =4 ??? Undefined function or variable "Ab".Error in ==> Gauss at 8[ap,p]=max(abs(Ab(k:n,k)));没有得到想要的结
果。
（2）第二题中在matlab中输入>> y=LU(A,b)得到y =4.0000 6.0000 -5.0000 -3.3571不是方程组的解。
（3）第三题中在用高斯赛德尔方法时在matlab中输入>> Gauss-Seidel(A,b,eps)结果程序报错??? Error using ==> Gauss Too many output arguments.得不到想要的结果。
4.2解决方案
（1）针对第一题中由于程序的第二行漏了一个分号导致输出了m和n的值，第8行中将Ab改为A问题就解决了。
（2）由于程序后面出现了矩阵y故输出的事矩阵y的值，但是我们要的事x的值，故只需要将y改成x，或者直接把y去掉就解决了问题。
（3）在function文件中命名不能出现“-”应该将其改为下划线“_”，所以将M文件名“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解决问题了。
五、心得体会
本次试验涉及到了用高斯列主元消去法，LU分解法，Jacobi迭代法以及Gauss-seidel迭代法等四种方法。需要对这些方法的原理都要掌握才能写出程序，由于理论知识的欠缺，我花了很大一部分时间在看懂实验的原理上，看懂了实验原理之后就开始根据原理编写程序，程序中还是出现了很多的低级错误导致调试很久才能运行。
通过这次试验使我深刻的体会到理论知识的重要性，没有理论知识的支撑是写不出程序来的。写程序时还会犯很多低级的错误，以后一定要加强理论知识的学习，减少编程时低级错误的产生。
牛顿迭代法文献综述
“牛顿迭代法”最新进展文献综述牛顿法是一种重要的迭代法，它是逐步线性化的方法的典型代表。牛顿迭代法又称为牛顿-拉夫逊方法，它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。另外该方法广泛用于计算机编程中。 介绍一下牛顿迭代法研究的前沿进展，1992年南京邮电学院基础课部的夏又生写的一篇题名一类代数方程组反问题的牛顿迭代法，对一类代数方程组反问题提出了一个可行的迭代解法。从算法上看,它是一种解正问题—迭代—解正问题迭代改善的求解过程。湖南师范大学的吴专保;徐大发表的题名堆浸工艺中浸润面的非线性问题牛顿迭代方法，为了研究堆浸工艺的机理,用牛顿迭代公式寻求浸润面的非线性方程的数值解,经过14次迭代的误差达到了,说明此算法收敛有效。浙江大学电机系的林友仰发表的牛顿迭代法在非线性电磁场解算中的限制对非线性电磁场解算中的限制做了分析，求解非线性方程组时迭代法是不可避免的。牛顿—拉斐森迭代法由于它的收敛速度快常被优先考虑。应用这个方法的主要问题是求雅可比矩阵。因为雅可比矩阵元素的计算非常费时。然而,本文要说明的是当利用以三角形为单元的有限元法求解非线性方程组时,应用牛顿法其雅可比矩阵容易求得,并且它保持了原系数的对称性和稀疏性,因而节省了时间。与此相反,若在差分法中应用牛顿迭代,并且按习惯用矩形网格进行剖分,则雅可比阵的计算很费时,而且不再保持原有对称性,这就使得存贮量和计算时间大为增加。南株洲工学院信息与计算科学系的吕勇;刘兴国发表的题名为牛顿迭代法加速收敛的一种修正格式，主要内容牛顿迭代法是求解非线性方程的一种重要的数值计算方法,在通常情况下,它具有至少平方收敛。本文利用文献[4]所建立的迭代格式xn+1=xn-αf(xfn)(x+n)f′(xn),对迭代格式中的参数α的讨论,实现了牛顿迭代法加速收敛的一种修正格式。
ICA使用牛顿迭代法对FastICA算法经行改进
ICA用牛顿迭代法改进的FastICA算法 ICA算法原理： 独立分量分析（ICA）的过程如下图所示：在信源()st中各分量相互独立的假设下，由观察xt通过结婚系统B把他们分离开来，使输出yt逼近st。 图1-ICA的一般过程 ICA算法的研究可分为基于信息论准则的迭代估计方法和基于统计学的代数方法两大类，从原理上来说，它们都是利用了源信号的独立性和非高斯性。基于信息论的方法研究中，各国学者从最大熵、最小互信息、最大似然和负熵最大化等角度提出了一系列估计算法。如FastICA算法, Infomax算法，最大似然估计算法等。基于统计学的方法主要有二阶累积量、四阶累积量等高阶累积量方法。本实验主要讨论FastICA算法。 1. 数据的预处理 一般情况下，所获得的数据都具有相关性，所以通常都要求对数据进行初步的白化或球化处理，因为白化处理可去除各观测信号之间的相关性，从而简化了后续独立分量的提取过程，而且，通常情况下，数据进行白化处理与不对数据进行白化处理相比，算法的收敛性较好。 若一零均值的随机向量 满足 ， 其中：I为单位矩阵，我们称这个向量为白化向量。白化的本质在于去相关，这同主分量分析的目标是一样的。在ICA中，对于为零均值的独立源信号 ， 有： ， 且协方差矩阵是单位阵cov( S ) = I，因此，源信号 S( t )是白色的。对观测信号X( t )，我们应该寻找一个线性变换，使X( t )投影到新的子空间后变成白化向量，即：
其中，W0为白化矩阵，Z为白化向量。 利用主分量分析，我们通过计算样本向量得到一个变换 其中U和 分别代表协方差矩阵XC的特征向量矩阵和特征值矩阵。可以证明，线性变换W0满足白化变换的要求。通过正交变换，可以保证 因此，协方差矩阵： 再将 代入 且令 有 由于线性变换A~连接的是两个白色随机矢量Z( t )和S( t )，可以得出A~ 一定是一个正交变换。如果把上式中的Z( t )看作新的观测信号，那么可以说，白化使原来的混合矩阵A简化成一个新的正交矩阵A~。证明也是简单的： 其实正交变换相当于对多维矢量所在的坐标系进行一个旋转。 在多维情况下，混合矩阵A是N*N 的，白化后新的混合矩阵A~ 由于是正交矩阵，其自由度降为N*(N-1)/2，所以说白化使得ICA问题的工作量几乎减少了一半。 白化这种常规的方法作为ICA的预处理可以有效地降低问题的复杂度，而且算法简单，用传统的PCA就可完成。用PCA对观测信号进行白化的预处理使得原来所求的解混合矩阵退化成一个正交阵，减少了ICA的工作量。此外，PCA本身具有降维功能，当观测信号的个数大于源信号个数时，经过白化可以自动将观测信号数目降到与源信号维数相同。
数值分析列主元消去法的实验报告
实验一 列主元消去法 【实验内容】 1．掌握列主元消去法的基本思路和迭代步骤 2．并能够利用列主元的高斯消去法解任意阶数的线性方程组； 3、从课后题中选一题进行验证，得出正确结果，交回实验报告与计算结果。 【实验方法与步骤】 1．列主元消去法基本思路 设有线性方程组Ax b =，设A 是可逆矩阵。列主元消去法的基本思想就是通过列主元的选取将初等行变换作用于方程组的增广矩阵[]|B A b =，将其中的A 变换成一个上三角矩阵，然后求解这个三角形方程组。 2．列主元高斯消去法算法描述 将方程组用增广矩阵[]()(1)|ij n n B A b a ?+==表示。 步骤1：消元过程，对1,2,,1k n =-L （1） 选主元，找{},1,,k i k k n ∈+L 使得 ,max k i k ik k i n a a ≤≤= （2） 如果,0k i k a =，则矩阵A 奇异，程序结束；否则执行（3）； （3） 如果k i k ≠，则交换第k 行与第k i 行对应元素位置，k kj i j a a ?， ,,1j k n =+L ； （4） 消元，对,,i k n =L ，计算/,ik ik kk l a a =对1,,1j k n =++L ，计算 .ij ij ik kj a a l a =- 步骤 2：回代过程： （1） 若0,nn a =则矩阵奇异，程序结束；否则执行（2）； （2） ,1/;n n n nn x a a +=对1,,2,1i n =-L ，计算 ,11/n i i n ij j ii j i x a a x a +=+??=- ??? ∑
[实验程序] #include
牛顿迭代法
牛顿迭代法 李保洋 数学科学学院信息与计算科学学号：060424067 指导老师：苏孟龙 摘要：牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法，即牛顿迭代法.迭代法是一种不断用变量的旧值递推新值的过程.跟迭代法相对应的是直接法或者称为一次解法，即一次性解决问题.迭代法又分为精确迭代和近似迭代.“牛顿迭代法”属于近似迭代法，本文主要讨论的是牛顿迭代法，方法本身的发现和演变和修正过程，避免二阶导数计算的Newton迭代法的一个改进，并与中国古代的算法，即盈不足术，与牛顿迭代算法的比较. 关键词：Newton迭代算法；近似求解；收敛阶；数值试验；中国古代数学； 九章算术；Duffing方程；非线性方程；收敛速度；渐进性 0 引言： 迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程，跟迭代法相对应的是直接法或者称为一次解法，即一次性解决问题.迭代法又分为精确迭代和近似迭代.“二分法”和“牛顿迭代法”属于近似迭代法. 迭代算法是用计算机解决问题的一种基本方法.它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令（或一定步骤）进行重复执行，在每次执行这组指令（或这些步骤）时，都从变量的原值推出它的一个新值.具体使用迭代法求根时应注意以下两种可能发生的情况： (1)如果方程无解，算法求出的近似根序列就不会收敛，迭代过程会变成死循环，因此在使用迭代算法前应先考察方程是否有解，并在程序中对迭代的次数给予限制. (2)方程虽然有解，但迭代公式选择不当，或迭代的初始近似根选择不合理，也会导致迭代失败. 所以利用迭代算法解决问题，需要做好以下三个方面的工作： 1、确定迭代变量.在可以用迭代算法解决的问题中，至少存在一个直接或间接地不断由旧值递推出新值的变量，这个变量就是迭代变量. 2、建立迭代关系式.所谓迭代关系式，指如何从变量的前一个值推出其下一个值的公式(或关系).迭代关系式的建立是解决迭代问题的关键，通常可以使用递推或倒推的方法来完成. 3、对迭代过程进行控制，在什么时候结束迭代过程？这是编写迭代程序必须考虑的问题.不能让迭代过程无休止地重复执行下去.迭代过程的控制通常可分为两种情况：一种是所需的迭代次数是个确定的值，可以计算出来;另一种是所需的迭代次数无法确定.对于前一种情况，可以构建一个固定次数的循环来实现对迭代过程的控制;对于后一种情况，需要进一步分析出用来结束迭代过程的条件. 1牛顿迭代法：
利用高斯列主元消去法求如下线性方程组的解
%利用高斯列主元消去法求如下线性方程组的解 clear all; A=[3 -2 1 -1;4 0 -1 2;0 0 2 3;0 0 0 5]; b=[8;-3;11;15]; function [X,XA] = UpGaussFun(A,b) %利用高斯列主元消去法求如下线性方程组的解 %A为一个n阶上三角非奇异矩阵 %b为线性方程组的阐述向量 %X为线性方程组AX=b的解 %XA为消元后的系数矩阵 N=size(A); n=N(1); index=0; for i=1:(n-1) me=max(abs(A(1:n,i)));%选列主元 for k=i:n if(abs(A(k,i))==me) index=k; break; end; end; end; temp=A(i,1:n); A(i,1:n)=A(index,1:n); A(index,1:n)=temp; bb=b(index); b(index)=b(i); b(i)=bb;%交换主行 for j=(i+1):n if(a(i,i)==0) disp('?????a???a0￡?'); return; end; l=A(j,i); m=A(i,i); A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; b(j)=b(j)-l*b(i)/m;
end; X=UpTriangleFun(A,b); XA=A; ----------------------------------------------------------------------------------------------------------------------------- % 函数定义 function [X,XA]= UpGaussFun(A,b) %利用高斯列主元消去法求如下线性方程组的解 %A为一个n阶上三角非奇异矩阵 %b为线性方程组的阐述向量 %X为线性方程组AX=b的解 %XA为消元后的系数矩阵 [N,M]=size(A); %N=sizes(A); n=N; index=0; for i=1:(n-1) me=max(abs(A(1:n,i))); %选列主元 for k=i:n if(abs(A(k,i))==me) index=k; break; end; end; temp=A(i,1:n); A(i,1:n)=A(index,1:n); A(index,1:n)=temp; bb=b(index); b(index)=b(i); b(i)=bb; %交换主行 for j=(i+1):n if(A(i,i)==0) disp('?????a???a0￡?'); return; end; l=A(j,i); m=A(i,i); A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; b(j)=b(j)-l*b(i)/m; end; end;
Newton迭代法实例
基于牛顿迭代法的圆形断面临界水深直接计算 学院：建筑工程学院学号：2111206052 姓名：王瑞峰 一、问题来源 圆形断面由于具有受力条件好、适应地形能力强、水力条件好等优点，已成为农田灌溉、城市给水排水等工程较常采用的断面形式。而临界水深的计算则是进行圆形断面水力计算的关键，但其计算较繁杂，要求解高次隐函数方程，且未知量包含在三角函数中，求解难度大。自20世纪90年代，对圆形断面临界水深的计算进行了大量研究，获得了较多成果。鉴此，本文应用牛顿迭代算法，得到一种较简洁且可提供高精度算法程序的近似计算公式。 二、数学模型 相应于断面单位能量最小值的水深称为临界水深，其计算公式为： 需满足的临界流方程为： 其中 式中，d为洞径；为临界水深对应的圆心角，rad；n为流速分布不均匀系数(不特殊说明时取1.0)；Q为流量，m3Is；g为重力加速度(通常取9.81 m/s2)；分别为临界流对应的过水断面面积和水面宽度。 无压流圆形断面的水力要素见图1 将式(1)、(3)、(4)代入式(2)得： 将式(5)整理即得临界水深的非线形方程： 由此可知．式(6)为临界水深h。的高次隐函数方程，且未知量包含在三角函数中。 即圆形断面临界水深的求解即为式(6)的求根问题。在现行工程实际中计算临界水深时均采用近似公式或试算法，所得结果精度不高且效率较低。 三、方法选择 牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。 解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点
附近展开成泰勒级数f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +… 取其线性部分，作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有f(x0)+f'(x0)(x－ x0)=f(x)=0 设f'(x0)≠0则其解为x1=x0－f(x0)/f'(x0) 这样，得到牛顿法的一个迭代序列：x(n+1)=x(n)－f(x(n))/f'(x(n))。 在对式(6)的求解方法中，应首选牛顿迭代法，因为牛顿迭代法可快速求解出其他方法求不出或难以求出的解。 引入无量纲参数k： 将式(7)代入式(6)得： 的一阶、二阶导函数分别为： 由牛顿迭代法可得： 式中，=0，1，2…为迭代次数；为的初值。 将式(8)、(9)代入式(10)，可得相应于式(6)临界水深对应中心角的牛顿迭代公式： 由式(11)迭代计算出临界水深对应的中心角后，代入式(1)即可得临界水深。 根据文献，为避免渡状水面有可能接触洞顶引起水流封顶现象。洞内水面线以上的空间不宜小于隧洞断面面积的15％，且高度不小于0.4m。可得临界水深对应的中心角的最大值一般不超过4.692，相应可得无量纲参数值的上限为0.5044。故取值范围为[O.000 0,0.504 4]。 查阅文献与的近似公式： 若将式(12)视为初值函数，代入式(11)进行一次迭代计算，不仅得到了直接计算的公式，且提高了计算结果的精度。 其中 将式(13)代入式(1)即得圆形断面临界水深。 计算实例： 某引水式电站输水隧洞为圆形断面，已知洞径d=3.0 m，试确定设计流量Q=8.0m3/s时的临界水深。 四、编程实现 本文采用Fortran软件求解，程序的代码如下：
各种迭代法编程
雅可比迭代法： function x=jacobi(a,b,p,delta,n) %a为n维非奇异矩阵；b为n维值向量 %p为初值；delta为误差界；n为给定的迭代最高次数 N=length（b）; for k=1:n for j=1：N x(j)=(b(j)-a(j,[1:j-1,j+1:N])*p([1:j-1,j+1:N]))/a(j,j); end err=abs(norm(x’-p)); p=x’; if(err function [x,k,err,p]=ddf(f,x0,tol,n) %ddl.m为用迭代法求非线性方程的解 %f为给定的迭代函数；x0为给定的初始值 %tol为给定的误差界；n为所允许的最大迭代次数 %k为迭代次数；x为不动点的近似值；err为误差 p(1)=x0; for k=2:n p(k)=feval(f,p(k-1)); k, err=abs(p(k)-p(k-1)) x=p(k); if(err 高斯列主元消元法解线性方程组 一、题目：用Gauss 列主元消去法解线性方程组Ax b =，其中， A=17.031 -0.615 -2.991 1.007 -1.006 0.000-1.000 34.211 -1.000 -2.100 0.300 -1.7000.000 0.500 13.000 -0.500 1.000 -1.5004.501 3.110 -3.907 -61.705 12.170 8.9990.101 -8.012 -0.017 -0.910 4.918 0.1001.000 2.000 3.000 4.500 5.000 21.803?? ? ? ? ? ? ? ? ??? 0.230 -52.322 54.000 240.236 29.304 -117.818b ?? ? ? ?= ? ? ? ? ??? T X=(0.907099 -1.961798 3.293738 -4.500708 3.029344 -5.255068) 二、原理及步骤分析 设 n n ij R a A ?∈=][)1(，n n R b b b b ∈=],,,[)1()2(2)1(1 。若约化主元素 ),,2,1(0)(n k a k kk =≠，则通过高斯消元法将方程b AX =约化为三角形方程组求解。 如果在消元过程中发现某个约化主元0) (=k kk a , 则第K 次消元就无法进行。此外，即 使所有约化主元全不为零，虽然可以完成方程组的求解，但也无法保证结果的可靠性，因为计算过程中存在舍入误差。 为减少计算过程中的舍入误差对解的影响，在每次消元前，应先选择绝对值尽可能大的元作为约元的主元，如果在子块的第一列中选取主元，则相应方法称为列主元消元法。相应过程为： （1）选主元：在子块的第一列中选择一个元) (k k i k a 使) (max k ik n i k k k i a a k ≤≤= 并将第k 行元与第k i 行元互换。 （2）消元计算：对k=1,2,……n-1依次计算 ()()()?? ?? ?????++=-=++=-=++==++n k k i b m b b n k k j i a m a a n k k i a a m k k ik k i k i k kj ik k ij k ij k kk k ik k ik ,,2,1,,2,1,,,2,1) ()()1() ()()1()() ()( （3）回代求解 用牛顿迭代法求近似根 ————————————————————————————————作者：————————————————————————————————日期： 第四题 题目：用Newton 法求方程在 74 28140x x -+= （0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）. 解：此题是用牛顿迭代法求解近似根的问题 1. Newton 迭代法的算法公式及应用条件： 设函数在有限区间[a,b]上二阶导数存在，且满足条件 ⅰ. ()()0f a f b <; ⅱ. ()''f x 在区间[a,b]上不变号； ⅲ. ()'0f x ≠; ⅳ. ()()'f c f c b a ≤-，其中c 是a,b 中使()()''min(,)f a f b 达到的一个. 则对任意初始近似值0[,]x a b ∈，由Newton 迭代过程 ()()() 1'k k k k k f x x x x f x +=Φ=-，k=0，1，2… 所生成的迭代序列{ k x }平方收敛于方程()0f x =在区间[a,b]上的唯一解а. 对本题： )9.1()9.1(0 )8(4233642)(0 )16(71127)(0 )9.1(,0)1.0(,1428)(3225333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f Θ 故以1.9为起点 ?? ???='-=+9.1)()(01x x f x f x x k k k k 2. 程序编写 #include 实验一 列主元消去法 【实验内容】1. 掌握列主元消去法的基本思路和迭代步骤 2. 并能够利用列主元的高斯消去法解任意阶数的线性方程组； 【实验方法与步骤】列主元消去法编写程序 1．列主元消去法基本思路 设有线性方程组Ax b =，设A 是可逆矩阵。列主元消去法的基本思想就是通过列主元的选取将初等行变换作用于方程组的增广矩阵[]|B A b =，将其中的A 变换成一个上三角矩阵，然后求解这个三角形方程组。 2．列主元高斯消去法算法描述 将方程组用增广矩阵[]()(1)|ij n n B A b a ?+==表示。 步骤1：消元过程，对1,2,,1k n =- （1） 选主元，找{},1,,k i k k n ∈+ 使得 ,max k i k ik k i n a a ≤≤= （2） 如果,0k i k a =，则矩阵A 奇异，程序结束；否则执行（3）； （3） 如果k i k ≠，则交换第k 行与第k i 行对应元素位置，k kj i j a a ?， ,,1j k n =+ ； （4） 消元，对,,i k n = ，计算/,ik ik kk l a a =对1,,1j k n =++ ，计算 .ij ij ik kj a a l a =- 步骤 2：回代过程： （1） 若0,nn a =则矩阵奇异，程序结束；否则执行（2）； （2） ,1/;n n n nn x a a +=对1,,2,1i n =- ，计算 ,11/n i i n ij j ii j i x a a x a +=+??=- ??? ∑ 习题3第一题程序如下 #include CENTRAL SOUTH UNIVERSITY 数值分析实验报告 非线性方程组的牛顿迭代法的应用 一、问题背景 非线性是实际问题中经常出现的，并且在科学与工程计算中的地位越来越重要，很多我们熟悉的线性模型都是在一定条件下由非线性问题简化的，为得到更符合实际的解答，往往需要直接研究非线性科学，它是21世纪科学技术发展的重要支柱，非线性问题的数学模型有无限维的如微分方程，也有有限维的。道遥咏计算机进行科学计算都要转化为非线性的单个方程或方程组的求解。从线性到非线性是一个质的变化，方程的性质有本质不同，求解方法也有很大差别。本文主要介绍的是非线性方程组的牛顿迭代法的数值解法。 二、数学模型 对于方程()0=x f ，如果()x f 湿陷性函数，则它的求根是容易的。牛顿法实质上是一种线性化方法，其基本思想是将线性方程()0=x f 逐步归结为某种线性方程来求解。 设已知方程()0=x f 有近似根k x （假定()0'≠k x f ），将函数()x f 在点k x 展开，有 ()()()()k k k x x x f x f x f -+≈'， 于是方程()0=x f 可近似地表示为 ()()()0'=-+k k k x x x f x f 这是个线性方程，记其根为1+k x ，则1+k x 的计算公式 ()() k k k k x f x f x x ' 1- =+， ,1,0=k 这就是牛顿法。 三、算法及流程 对于非线性方程 ()()()???? ????????=n n n n x L x x f M x L x x f x L x x f f ,,,,,,,,,2 1212211 在()k x 处按照多元函数的泰勒展开，并取线性项得到 非线性回归预测法 前面所研究的回归模型，我们假定自变量与因变量之间的关系是线性的，但社会经济现象是极其复杂的，有时各因素之间的关系不一定是线性的，而可能存在某种非线性关系，这时，就必须建立非线性回归模型。 一、非线性回归模型的概念及其分类 非线性回归模型，是指用于经济预测的模型是曲线型的。常见的非线性回归模型有下列几种： (1)双曲线模型： i i i x y εββ++=1 2 1 (3-59) (2)二次曲线模型： i i i i x x y εβββ+++=2321 (3-60) (3)对数模型： i i i x y εββ++=ln 21 (3-61) (4)三角函数模型： i i i x y εββ++=sin 21 (3-62) (5)指数模型： i x i i ab y ε+= (3-63) i i i x x i e y εβββ+++=221110 (3-64) (6)幂函数模型： i b i i ax y ε+= (3-65) (7)罗吉斯曲线： i x x i i i e e y εββββ++=++1101101 (3-66) (8)修正指数增长曲线： i x i i br a y ε++= (3-67) 根据非线性回归模型线性化的不同性质，上述模型一般可细分成三种类型。 第一类：直接换元型。 这类非线性回归模型通过简单的变量换元可直接化为线性回归模型，如：(3-59)、(3-60)、(3-61)、(3-62)式。由于这类模型的因变量没有变形，所以可以直接采用最小平方法估计回归系数并进行检验和预测。 第二类：间接代换型。 这类非线性回归模型经常通过对数变形的代换间接地化为线性回归模型，如：(3-63)、(3-64)、(3-65)式。由于这类模型在对数变形代换过程中改变了因变量的形态，使得变形后模型的最小平方估计失去了原模型的残差平方和为最小的意义，从而估计不到原模型的最佳回归系数，造成回归模型与原数列之间的较大偏差。 第三类：非线性型。 实验三 高斯列主元消去法 一、实验目的： 1、掌握高斯消去法的基本思路和迭代步骤。 2、 培养编程与上机调试能力。 二、高斯列主元消去法的基本思路与计算步骤： 设有方程组Ax b =，设A 是可逆矩阵。高斯消去法的基本思想就是僵局真的初等行变换作用于方程组的增广矩阵[]B A b = ，将其中的A 变换成一个上三角矩阵，然后求解这个三角形方程组。 列主元高斯消去法计算步骤： 将方程组用增广矩阵[]()(1)ij n n B A b a ?+== 表示。 步骤1：消元过程，对1,2,,1k n =- （1） 选主元，找{},1,,k i k k n ∈+ 使得 ,max k i k ik k i n a a ≤≤= （2） 如果 ,0k i k a =，则矩阵A 奇异，程序结束；否则执行（3）。 （3） 如果k i k ≠，则交换第k 行与第k i 行对应元素位置，k kj i j a a ?，,,1j k n =+ 。 （4） 消元，对,,i k n = ，计算/,ik ik kk l a a =对1,,1j k n =++ ，计算 . ij ij ik kj a a l a =- 步骤 2：回代过程： （1） 若0,nn a =则矩阵奇异，程序结束；否则执行（2）。 （2） ,1/;n n n nn x a a +=对1,,2,1i n =- ，计算,11/n i i n ij j ii j i x a a x a +=+??=- ???∑ 三：程序流程图 四：程序清单： function X=uptrbk(A,b) % A是一个n阶矩阵。 % b是一个n维向量。 % X是线性方程组AX=b的解。 [N N]=size(A); X=zeros(1,N+1); Aug=[A b]; for p=1:N-1 [Y,j]=max(abs(Aug(p:N,p)));%返回向量的最大值存入y，最大值的序号存入j。 C=Aug(p,:); Aug(p,:)=Aug(j+p-1,:); Aug(j+p-1,:)=C; if Aug(p,p)==0 'A是奇异阵，方程无惟一解' break end for k=p+1:N m=Aug(k,p)/Aug(p,p); Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1); end end % 这里用到程序函数backsub来进行回代。 X=backsub(Aug(1:N,1:N),Aug(1:N,N+1)); function X=backsub(A,b) % A是一个n阶上三角非奇异阵。 % b是一个n维向量。 % X是线性方程组AX=b的解。 n=length(b);%取b向量的个数。 X=zeros(n,1); X(n)=b(n)/A(n,n); for k=n-1:-1:1 X(k)=(b(k)-A(k,k+1:n)*X(k+1:n))/A(k,k); End 五、测试数据与结果： 测试数据：（第8章习题三第2题）求解线性方程组： 解：建立一个主程序gs.m clc clear A=[1,2,3;5,4,10;3,-0.1,1]; b=[1;0;2]; 1.填空 1) 计算 f=(2-1)6 , 取2＝1.4 , 利用下列算式，那个得到的结果最好？答：C (A) 6121 )(-, (B) (3-22)2, (C) 32231)(+, (D) 99-702 2) 称序列{x n }是p 阶收敛的条件为c x x x x p n n n =--+∞→** lim 1 3) 在等式∑==n k k k n x f a x x x f 010)(],,,[ 中, 系数a k 与函数f (x ) 无 关。 （限填“有”或“无”） 4) 设P k (x k ,y k ) , k =1,2,…,5 为函数y =x 2-3x +1上的5个互异的点，过P 1,…,P 5且次数不超过4次的插值多项式是 x 2-3x +1 。 5) 设f (x )∈C [a ,b ]， f (x )的最佳一致逼近多项式是__一定___存在的。 6) 求解微分方程数值解的E ul e r 法的绝对稳定区间是(-2,0) 。 7) n 个节点的插值型求积公式的代数精度不会超过2n －1次。 8) 高次插值容易产生________龙格（R u n g e ）现象。 9) R n 上的两个范数||x||p , ||x||q 等价指的是_?C,D ∈R,_C_||x||q _≤||x||p ≤D ||x||q _； R n 上的两个范数_一定__是等价的。（选 填“一定”或“不一定”）。 2.曲线151.03+-=x x y 与89.14.22-=x y 在点（1.6,1）附近相切，试用牛顿迭代法求切点横坐标的近似值1+k x ，使5110-+≤-k k x x 。 解 两曲线的导数分别为51.032-='x y 和x y 8.4='，两曲线相切，导数相等，故有 051.08.432=--x x 令51.08.43)(2--=x x x f ，则f(1)<0,f(2)>0，故区间[1,2]是f(x)=0的有根区间，又当]2,1[∈x 时，08.46)(>-='x x f ，因此f(x)=0在[1,2]上有惟一实根x*，对f(x)应用牛顿迭代法，得计算公式 ,2,1,0,8 .4651.08.4321=----=+k x x x x x k k k k k 由于06)(>=''x f ，故取20=x 迭代计算一定收敛，计算结果如表7-6所示。 表7-6 k k x k k x 0 2.0 3 1.706815287 1 2.293055556 4 1.700025611 2 1.817783592 5 1.7 继续计算仍得7.16=x ，故7.1*=x 。 注 本题也可令89.14.2151.02 3-=+-x x x ，解得切点横坐标满足方程089.2514.2)(23=+--=x x x x f ，用有重根时的牛顿迭代法（7.15）式计算，此时m=2，仍取x0=2，经四步可得x*=1.7。 clear;clc; %Gauss消去法解线性方程组 A=[3 -5 6 4 -2 -3 8; 1 1 -9 15 1 -9 2; 2 -1 7 5 -1 6 11; -1 1 3 2 7 -1 -2; 4 3 1 -7 2 1 1; 2 9 -8 11 -1 -4 -1; 7 2 -1 2 7 -1 9];%系数矩阵 b=[11 2 29 9 5 8 25]';%n维向量 y=inv(A)*b %matlab的计算结果 n=length(b);%方程个数n x=zeros(n,1);%未知向量 %-------------消去----------- for k=1:n-1 % if A(k,k)==0; % error('Error'); % end for i=k+1:n % A(i,k)=A(i,k)/A(k,k); Aik=A(i,k)/A(k,k) for j=k:n A(i,j)=A(i,j)-Aik*A(k,j); end A b(i)=b(i)-Aik*b(k) end end %-------------回代----------- x(n)=b(n)/A(n,n) for k=n-1:-1:1 S=b(k); for j=k+1:n S=S-A(k,j)*x(j); end x(k)=S/A(k,k) end x %程序的计算结果 error=abs(x-ones(n,1))%误差 clear;clc; %列主元Gauss校区法解线性方程组 A=[3 -5 6 4 -2 -3 8; 1 1 -9 15 1 -9 2; 2 -1 7 5 -1 6 11; -1 1 3 2 7 -1 -2; 4 3 1 -7 2 1 1; 2 9 -8 11 -1 -4 -1; 7 2 -1 2 7 -1 9];%系数矩阵 b=[11 2 29 9 5 8 25]';%n维向量 y=inv(A)*b %matlab的计算结果 n=length(b);%方程个数n x=zeros(n,1);%未知向量 %-------------消去----------- for k=1:n-1 Auk=A(k:n,k); [m,u]=max(abs(Auk)); u=u+k-1 %u为最大元所在的列 %------交换最大的行和当前行的值------- for j=k:n temp=A(u,j);A(u,j)=A(k,j);A(k,j)=temp; end temp=b(k);b(k)=b(u);b(u)=temp; % if A(k,k)==0; % error('Error'); % end for i=k+1:n % A(i,k)=A(i,k)/A(k,k); Aik=A(i,k)/A(k,k) for j=k:n A(i,j)=A(i,j)-Aik*A(k,j); end A b(i)=b(i)-Aik*b(k) end end %-------------回代----------- x(n)=b(n)/A(n,n) for k=n-1:-1:1 S=b(k); for j=k+1:n S=S-A(k,j)*x(j); 高斯消元法1.程序: clear format rat A=input('输入增广矩阵A=') [m,n]=size(A); for i=1:(m-1) numb=int2str(i); disp(['第',numb,'次消元后的增广矩阵']) for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); end A end %回代过程 disp('回代求解') x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); end x 2.运行结果： 高斯选列主元消元法1.程序: clear format rat A=input('输入增广矩阵A=') [m,n]=size(A); for i=1:(m-1) numb=int2str(i); disp(['第',numb,'次选列主元后的增广矩阵']) temp=max(abs(A(i:m,i))); [a,b]=find(abs(A(i:m,i))==temp); tempo=A(a(1)+i-1,:); A(a(1)+i-1,:)=A(i,:); A(i,:)=tempo disp(['第',numb,'次消元后的增广矩阵']) for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); end A end %回代过程 disp('回代求解') x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); end x 2.运行结果： 《Gauss列主元消去法》实验报告 实验名称：Gauss列主元消去法程序设计？？？成绩：_________ 专业班级：数学与应用数学1202班？姓名：王晓阳？??学号： 实？验？日？期：？2014?年11月10日 实验报告日期：？2014年？11月10日 一.实验目的 1. 学习Gauss消去法的基本思路和迭代步骤. 2. 学会运用matlab编写高斯消去法和列主元消去法程序，求解线性方程组. 3. 当绝对值较小时，采用高斯列主元消去法? 4. 培养编程与上机调试能力. 二、实验内容 用消去法解线性方程组的基本思想是用逐次消去未知数的方法把原线性方程组Ax二b 化为与其等价的三角形线性方程组，而求解三角形线性方程组可用回代的方法求解 1. 求解一般线性方程组的高斯消去法? (1) 消元过程： 设a kk k-0 ,第i个方程减去第k个方程的m ik Tk k倍，("k 1^1, n)，得到 A k1x=b k1. 经过n-1次消元，可把方程组A1^b1化为上三角方程组A n x=b n. ⑵回代过程: 以解如下线性方程组为例测试结果 2. 列主元消去法 由高斯消去法可知，在消元过程中可能出现a kk k =0的情况，这是消去法将无法进行, 即使主元素a kk k-0但很小时，用其作除数，会导致其他元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠.这时就需要选取主元素，假定线性方程组的系数矩阵A是菲奇异的. (1)消元过程： 对于k =1,2,川,n -1，进行如下步骤: 1) 按列选主元，记 2) 交换增广阵A的p,k两行的元素 A(k,j)=A(p,j) ( j=k,…,n +1) 3) 交换常数项b的p,k两行的元素。 b(k)=b(p) 4) 计算消元 (2) 回代过程 (3) 以解如下线性方程组为例测试结果 三、实验环境 MATLAB R2014a 四、实验步骤高斯列主元消元法解线性方程组
用牛顿迭代法求近似根
列主元消去法
非线性方程组的牛顿迭代法的应用
非线性回归预测法——高斯牛顿法(詹学朋)
实验三高斯列主元消去法
数值分析练习第五套
高斯法和列主元高斯消去法解线性方程组(MATLAB版)
数值分析编程及运行结果(高斯顺序消元法)
Gauss列主元消去法程序设计
- 列主元消去法解方程组实验报告
- Gauss列主元消去法程序设计
- 高斯列主元消去法
- 高斯列主元消元法解线性方程组
- Gauss列主元消去法实验报告
- 列主元高斯消去法实验报告
- Gauss列主元消去法程序设计
- 计算方法-实验三列主元高斯消去法
- 高斯列主元消去法
- 数值分析列主元消去法的实验报告
- LU分解高斯消元列主元高斯消元matlab代码
- 数值分析列主元消去法的实验报告
- 高斯列主元消去法解线性方程组VC++课程设计报告
- 高斯列主元消元法解线性方程组
- matlab高斯列主元消去法
- 选列主元的高斯消去法实验报告2
- 实验二：列主元消元法实验报告
- Gauss列主元消去法程序设计说明书
- 列主元高斯约当消去法实验讲解
- Gauss列主元消去法程序设计.doc