(完整版)第八节雅可比迭代法与高斯-塞德尔迭代法综述
第八节 雅可比迭代法与高斯—塞德尔迭代法
一 雅可比迭代法
设线性方程组
b Ax = (1) 的系数矩阵A 可逆且主对角元素nn a ,...,a ,a 22
11均不为零,令
()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)迭代法(公式),用向量的分量来表示,(4)为
?
??
[]
,...
,,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=为初始向量.
由此看出,雅可比迭代法公式简单,每迭代一次只需计算一次矩阵和向量的乘法.在电算时需
要两组存储单元,以存放()
k x 及()
1+k x . 例1 例1 用雅可比迭代法求解下列方程组
???
??=+--=-+-=--2
45382102
7210321321321.x x x .x x x .x x x
解 将方程组按雅可比方法写成
??
??
???++=++=++=8402020830201072
020*******
2321.x .x .x .x .x .x .x .x .x
取初始值
()()()()
()()T T ,,,x ,x ,x x 0000302010==按迭代公式
()()
()()()
()()()()???????++=++=++=+++840202083020107202010211331
123211.x .x .x .x .x .x .x .x .x k k k k k k k k k
进行迭代，其计算结果如表1所示
表1
二 高斯—塞德尔迭代法
由雅可比迭代公式可知,在迭代的每一步计算过程中是用()
k x
的全部分量来计算()
1+k x
的所
有分量,显然在计算第i 个分量()
1+k i x 时,已经计算出的最新分量()
()
11
11
+-+k i k x ,...,x 没有被利
用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第
1+k 次近似()
1+k x
的分量
()
1+k j
x 加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）
迭代法.
把矩阵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==∑∑--=-=+=++Λ (9)
由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出
()
1+k i
x 后()
k i
x 不再使用,所以用()
1+k i x 冲掉()
k i
x ,以便存放近似解.
例2 例2 用高斯——塞德尔迭代法求解例1.
解 取初始值
()()()()
()()T
T
,,,x ,x ,x x 0000302010==，按迭代公式
()()
()()()
()()()()
?????++=++=++=++++++840202083020107202010121113
311
123211.x .x .x .x .x .x .x .x .x k k k k k k k k k
进行迭代，其计算结果如下表2
从此例看出,高斯—塞德尔迭代法比雅可比迭代法收敛快(达到同样的精度所需迭代次数少),但这个结论,在一定条件下才是对的,甚至有这样的方程组,雅可比方法收敛，而高斯—塞德尔迭代法却是发散的.
三 迭代收敛的充分条件
定理1 在下列任一条件下,雅可比迭代法(5)收敛.
①
1
11
<∑
=≠=∞
n
i
j j ii
ij i
a a max B ;
②
1
11
1
<∑
=≠=n
i
j i ii
ij j
a a max B ;
③ 1
11
<∑
=-≠=∞
-n
j
i i jj
ij j
T
a a max A
D I
定理2 设
21B B ,分别为雅可比迭代矩阵与高斯—塞德尔迭代矩阵,则
∞
∞
≤1
2
B B (10)
从而，当
1
11
<∑
=≠=∞
n
i
j j ii
ij i
a a max B
时，高斯—塞德尔迭代法(8)收敛. 证明 由
21B B ,的定义,它们可表示成
()U L D B +=-11
()()U D L D I U L D B 1111
2-----=-=
用e 表示n 维向量()T
,...,,e 111=,则有不等式
e
B e B ∞≤11
U
D L D B 111--+=
这里,记号｜·｜表示其中矩阵的元素都取绝对值,而不等式是对相应元素来考虑的,于是
()()()Ie B L D I e
L D B e U D ∞
------≤-=1
11111
容易验证
()
1
1
==--n
n
L D L D
所以,L D I 1
--及
L D I 1
--可逆，且
(
)
()()
1
11
111
111
1-----------=++≤+++=-L D I L
D ...L D I L D ...L D I L
D I n n
()
I L D I ≥---1
1
从而有
()
()((){}e I B L D I L D I e
U D L
D I e B ∞
----------≤?-≤11
1
111
121
{()()}e
B e
L D I I B I ∞--∞
≤-?--=11
1
11
因此必有
∞
∞≤1
2
B B
因为已知
11
<∞
B 所以12<∞
B .即高斯—塞德尔迭代法收敛.
若矩阵
A 为对称，我们有
定理3 若矩阵A 正定,则高斯—塞德尔迭代法收敛.
证明 把实正定对称矩阵A 分解为
T L L D A --=
()T
L U =,则D 为正定的,迭代矩阵
()T L L D B 1
2--=
设λ是
2B 的任一特征值,x 为相应的特征向量,则
()
()x x L L D T λ=--1
以L D -左乘上式两端,并由T
L L D A --=有
()Ax x L T λλ=-1
用向量x 的共轭转置左乘上式两端,得
()Ax x x L x
T
T
T
--=-λλ1 (11)
求上式左右两端的共轭转置,得
Ax x x L x T T -
-
-
-=??? ??-λλ1
以λ-
-1和λ-1分别乘以上二式然后相加,得
()()Ax
x x L L x T T T -
-
---??? ??-+=+??? ??--λλλλλλ211 由T
L L D A --=,得
()()Ax
x x A D x T T -----??? ??-+=-??? ??--λλλλλλ211
即
(
)
Ax x x L x T
T
--
-=-λλλ2
2
11 (12)
因为A 和D 都是正定的,且x 不是零向量,所以由(11)式得1≠λ
,而由(12)式得
012
>-λ, 即1<λ,从而()12
定义1 设)n
n ij a A ?=为n 阶矩阵.
① ①如果
n
,...,i ,
a a n
i
j i j ij ii 21=∑>≠= (13)
即A 的每一行对角元素的绝对值都严格大于同行其他元素绝对值之和，则称A 为严格对角优势矩阵.
② ②如果
n
,...,i ,
a a n
i
j i j ij ii 21=∑≥≠=
且至少有一个不等式严格成立,则称A 为弱对角优势矩阵.
例如?????????
?-310131012是严格对角优势矩阵，??????????--310121011是弱对角优势矩阵. 定义2 设()n
n ij a A ?=是n 阶矩阵，如果经过行的互换及相应列的互换可化为???
??
?221211
0A A A ,
即存在n 阶排列矩阵P,使
???
???=2212110A A A AP P T
其中2211A ,A 为方阵，则称A 是可约的,否则称A 为不可约的.
A 是可约矩阵,意味着b Ax =可经过若干次行列重排,化为两个低阶方程组,事实上,
b Ax =可化为 ()b P x P AP P T
T T =,记
()()()()?
?????==??????==2121d d d b P ,y y y x P T
T
于是，求解b Ax =化为求解
()()()()()
?????=+=+22221212111d y A d y A y A
可以证明,如果A 为严格对角优势矩阵或为不可约弱对角优势矩阵,则A 是非奇异的.
定理4 如果A 为严格对角优势矩阵或为不可约弱对角优势矩阵,则对任意()
0x ,雅可比迭代法(4)与高斯—塞德尔迭代法(8)均为收敛的.
证明 下面我们以A 为不可约弱对角优势矩阵为例，证明雅可比迭代法收敛，其他证明留给读者.
要证明雅可比迭代法收敛，只要证()11
用反证法,设矩阵
1B 有某个特征值μ,使得1≥μ,则()01=-B I det μ,由于A 不可约,
且具有弱对角优势，所以1
-D 存在，且 ()()D A D D A D I I B I -+=--=---μμμ111
从而
()0=-+D A D det
μ
另一方面，矩阵()D A D -+μ与矩阵A 的非零元素的位置是完全相同的，所以()D A D -+μ也是不可约的,又由于1≥μ,且A 弱对角优势，所以
n
,...,i ,
a a a n
i
j i j ij ii ii 21=∑≥≥≠=μ
并且至少有一个i 使不等号严格成立.因此,矩阵
()D A D -+μ弱对角优势，故
()D A D -+μ为不可约弱对角优势矩阵.从而
()0≠-+D A D det μ
矛盾，故
1B 的特征值不能大于等于1，定理得证.
高斯-赛德尔迭代法matlab程序
disp('划分为M*M个正方形') M=5 %每行的方格数，改变M可以方便地改变剖分的点数 u=zeros(M+1);%得到一个（M+1）*（M+1）的矩阵 disp('对每个剖分点赋初值，因为迭代次数很高，所以如何赋初值并不重要，故采用对列线性赋值。') disp('对边界内的点赋初值并使用边界条件对边界赋值：') for j=1:M-1 for i=1:M-1 u(i+1,j+1)=100*sin(pi/M*j)/M*(M-i);%对矩阵（即每个刨分点）赋初值 end end for i=1:M+1 u(1,i)=100*sin(pi*(i-1)/M);%使用边界条件对边界赋值 u(1,M+1)=0; end u tic %获取运行时间的起点 disp('迭代次数为N') N=6 %迭代次数，改变N可以方便地改变迭代次数 disp('n为当前迭代次数，u为当前值，结果如下：') for n=1:N for p=2:M i=M+2-p; for j=2:M u(i,j)=0.25*(u(i,j-1)+u(i+1,j)+u(i-1,j)+u(i,j+1));%赛德尔迭代法 end end n %输出n u %输出u end disp('所用的时间：') t=toc %获取算法运行需要的时间 [x,y]=meshgrid(0:1/M:1,0:1/M:1); z=u(1,:); for a=2:M+1 z=[z;u(a,:)];%获取最终迭代的结果，幅值给z，z的值代表该点的点位值 end mesh(x,y,z)%绘制三维视图以便清楚地显示结果 mesh(x,y,z,'FaceColor','white','EdgeColor','black') %绘制三维视图以便清楚地显示结果
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本身具有降维功能，当观测信号的个数大于源信号个数时，经过白化可以自动将观测信号数目降到与源信号维数相同。
雅可比迭代法
2013-2014（1）专业课程实践论文 题目：雅可比迭代法
一、算法理论 设有方程组),...,2,1(1 n i b x a i j n j ij ==∑= 记作,b Ax = （1） A 为非奇异阵且),,...,2,1(0n i a ij =≠将A 分裂为U L D A --=，其中 D =????????????????nn a a a 22 11，L =-??? ????? ???? ????-00001,21323121n n n n a a a a a a U =-?? ? ?? ? ? ? ????????-0000,122311312n n n n a a a a a a 将式（1）第)....2,1(n i i =个方程用ii a 去除再移项，得到等价方程组 (),,...,2,111n i x a b a x n i j j j ij i ii i =??? ? ? ?? -=∑≠= （2） 简记作 ,0f x B x += 其中 ().,111 0b D f U L D A D I B ---=+=-= 对方程组（2）应用迭代法，得到解式（1）的雅可比迭代公式 () () ()()()()()????????? ?? ? ??- ==∑≠=+,1,...,11002010n i j i k j ij i ii k i t n x a b a x x x x x ， 初始向量 （3）
其中()()()()()T k n k k k x x x x ,,...,21=为第k 次迭代向量。设()k x 已经算出，由式（3）可计算下一次迭代向量()(),,...,2,1,...;2,1,01n i k x k ==+ 显然迭代公式（3）的矩阵形式为 ()()()()???+=+,010f x B x x k k ，初始向量 其中0B 称为雅可比方法迭代矩阵。
MATLAB样例之雅克比迭代法
要求： 下面分别使用雅克比迭代法和高斯-赛德尔迭代法求一个方程组的近似解用的线性方程组是按实验要求给的： 7*x1+x2+2*x3=10 x1+8*x2+2*x3=8 2*x1+2*x2+9*x3=6 雅克比迭代法的matlab代码：（老师写的） A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(any(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); while 1 x1=B*x0+f K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf)
A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(all(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); x00=x0; while 1 x11=B*x0+f; x00(1,1)=x11(1,1); x12=B*x00+f; x00(2,1)=x12(2,1); x13=B*x00+f; x00(3,1)=x13(3,1); x1=x00 K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf)
各种迭代法编程
雅可比迭代法： 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一Seidel method)是求解电力系统潮流的方法。潮流计算高斯——赛德尔迭代法又分导纳矩阵迭代法和阻抗矩阵迭代法两种。前者是以节点导纳矩阵为基础建立的赛德尔迭代格式;后者是以节点阻扰矩阵为基础建立的赛德尔迭代格式。高斯——赛德尔迭代法这是数学上求解线性或非线性方程组的一种常用的迭代方法。 本实验通过对电力网数学模型形成的计算机程序的编制与调试，获得形成电力网数学模型：高斯---赛德尔法的计算机程序，使数学模型能够由计算机自行形成，即根据已知的电力网的接线图及各支路参数由计算程序运行形成该电力网的节点导纳矩阵和各节点电压、功率。通过实验教学加深学生对高斯---赛德尔法概念的理解，学会运用数学知识建立电力系统的数学模型，掌握数学模型的形成过程及其特点，熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。 高斯---赛德尔法潮流计算框图 ［1］系统节点的分类 根据给定的控制变量和状态变量的不同分类如下 ①P、Q节点（负荷节点），给定Pi、Qi求Vi、Si，所求数量最多； ②负荷节点，变电站节点（联络节点、浮游节点），给定P Gi、Q Gi的发电机 节点，给定Q Gi的无功电源节点； ③PV节点（调节节点、电压控制节点），给定P i、Q i求Q n、S n，所求数量 少，可以无有功储备的发电机节点和可调节的无功电源节点； ④平衡节点（松弛节点、参考节点（基准相角）、S节点、VS节点、缓冲节 点），给定V i，δi=0，求P n、Q n(V s、δs、P s、Q s)。 ［2］潮流计算的数学模型 1)线性的节点电压方程YV=I 根据S=V错误!未找到引用源。可得非线性的节点电压方程(错误!未找到引用源。为I的共轭) YV=I=错误!未找到引用源。=错误!未找到引用源。 非线性回归预测法 前面所研究的回归模型，我们假定自变量与因变量之间的关系是线性的，但社会经济现象是极其复杂的，有时各因素之间的关系不一定是线性的，而可能存在某种非线性关系，这时，就必须建立非线性回归模型。 一、非线性回归模型的概念及其分类 非线性回归模型，是指用于经济预测的模型是曲线型的。常见的非线性回归模型有下列几种： (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)式。由于这类模型在对数变形代换过程中改变了因变量的形态，使得变形后模型的最小平方估计失去了原模型的残差平方和为最小的意义，从而估计不到原模型的最佳回归系数，造成回归模型与原数列之间的较大偏差。 第三类：非线性型。 数值分析实验五 班级： 10信计二班 学号：59 姓名：王志桃 分数： 一．实验名称 高斯-赛德尔迭代法解线性方程组 二．实验目的 1. 学会利用高斯赛德尔方法解线性方程组 2. 明白迭代法的原理 3. 对于大型稀疏矩阵方程组适用于迭代法比较简单 三．实验内容 利用Gauss-Seidel 迭代法求解下列方程组 ?????=++=-+=+-36123633111420238321 321321x x x x x x x x x ， 其中取→=0)0(x 。 四、算法描述 由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值，若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量)1(+k i x 时，用最新分量)1(1+k x ，???+)1(2k x )1(1-+k i x 代替旧分量)(1k x ，???)(2k x )(1-k i x ，就得到所谓解方程组的Gauss-Seidel 迭代法。 其迭代格式为 T n x x x x ) ()0()0(2)0(1)0(，，，???= (初始向量)， )(1111 1) ()1()1(∑∑-=-+=++--=i j i i j k j ij k j ij i ii i i x a x a b a x )210i 210(n k ???=???=，，，；，，， 或者写为 ?? ???--=???=???==?+=∑∑-=-+=+++)(1)210i 210(1111)()1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，， 五、 编码 #include 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。 高斯消元法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.运行结果： 分别运用高斯赛德尔迭代法和超松弛迭代法解线性方程组：????? ??-=????? ??????? ??--243024410143034321x x x 。 1． 高斯赛德尔迭代法 编程思路：高斯赛德尔迭代法实在雅克比迭代法的基础上进行优化得到的，即在进行迭代时，将已经算得的第k+1步的迭代值代入第k+1步后边的变量的计算当中去，从而加快了迭代速度。 程序代码： function varargout=Gauss_Seidelli(varargin) A=[4 3 0;3 4 -1;0 -1 4]; b=[24 30 -24]'; x0=[0;0;0]; x=Gauss_Seidel(A,b,x0) function x=Gauss_Seidel(A,b,x0) n=100;%最大迭代次数 ee=0.0001;%精度 n1=length(b); for i=1:n x1=x0; for j=1:n1 s=0; for k=1:n1 if k~=j s=s+A(j,k)*x0(k); end end x0(j)=(b(j)-s)/A(j,j); end if norm(x1-x0) 2． 超松弛迭代法 该方法是在高斯赛德尔迭代法的基础上将前一步的结果)(k i x 和)1( k i x 进行适当的线性组合以加速收敛，松弛因子ω的选择是关键，当1<ω<2时，即为超松弛迭代法。 程序代码： function varargout=SORli(varargin) clc A=[4 3 0;3 4 -1;0 -1 4]; b=[24;30;-24]; x0=[0;0;0];w=1.3; x=SOR(A,b,x0,w); for i=1:3 fprintf('%4.2f ',x(i)); end fprintf('\n'); function x=SOR(A,b,x0,w) %AX=b %x0初始点 %w 为 松弛因子 n=100;%最大迭代次数 ee=0.0001;%精度 n1=length(b); for i=1:n x1=x0; for j=1:n1 s=0; for k=1:n1 if k~=j s=s+A(j,k)*x0(k); end end x0(j)=(1-w)*x0(j)+w*(b(j)-s)/A(j,j); end if norm(x1-x0) 例1、已知函数表 求() f x的Lagrange二次插值多项式和Newton二次插值多项式。 解： （1）由题可知 插值基函数分别为 故所求二次拉格朗日插值多项式为 （2）一阶均差、二阶均差分别为 均差表为 故所求Newton 二次插值多项式为 例2、 设2()32f x x x =++，[0,1]x ∈，试求()f x 在[0, 1]上关于()1x ρ=，{}span 1,x Φ=的 最佳平方逼近多项式。 解： 若{}span 1,x Φ=，则0()1x ?=，1()x x ?=，且()1x ρ=，这样，有 所以，法方程为 011231261192 34a a ??????????=?????????? ?????????? ，经过消元得012311 62110123a a ??? ???????=???????????????????? 再回代解该方程，得到14a =，011 6 a = 故，所求最佳平方逼近多项式为* 111 ()46 S x x = + 例3、 设()x f x e =，[0,1]x ∈，试求()f x 在[0, 1]上关于()1x ρ=，{}span 1,x Φ=的最佳平 方逼近多项式。 解： 若{}span 1,x Φ=，则0()1x ?=，1()x x ?=，这样，有 所以，法方程为 解法方程，得到00.8732a =，1 1.6902a =， 故，所求最佳平方逼近多项式为 例4、 用4n =的复合梯形和复合辛普森公式计算积分1?。 解： （1）用4n =的复合梯形公式 由于 2h =，()f x =，()121,2,3k x k k =+=，所以，有 （2）用4n =的复合辛普森公式 由于2h =，()f x =，()121,2,3k x k k =+=，()1 2 220,1,2,3k x k k +=+=，所以，有 例5、 用列主元消去法求解下列线性方程组的解。 解：先消元 再回代，得到33x =，22x =，11x = 实验五 矩阵的lu分解法，雅可比迭代法 班级： 学号： 姓名： 实验五 矩阵的LU 分解法，雅可比迭代 一、目的与要求： 熟悉求解线性方程组的有关理论和方法； 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序； 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。 二、实验内容： 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解 各种方法的优缺点。 三、程序与实例 列主元高斯消去法 算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+?表示 1) 消元过程 对k=1,2,…,n-1 ①选主元，找{}n ,,1k ,k i k +∈使得 k ,i k a = ik a n i k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。 ③如果k i k ≠，则交换第k 行与第k i 行对应元素位置， j i kj k a a ? j=k,┅,n+1 ④消元，对i=k+1, ┅,n 计算 kk ik ik a a l /= 对j=l+1, ┅,n+1计算 kj ik ij ij a l a a -= 2) 回代过程 ①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。 ②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算 ii n i j j ij n i i a x a a x /11,??? ? ? ?- =∑+=+ 程序与实例 程序设计如下： #include 一、 实验目的与要求 对于线性方程组?????=++=++=++69228281027321 321321x x x x x x x x x 1. 用高斯-赛德尔迭代法求此方程组的近似解（终止迭代过程的最大允许迭代次数N ，近似解的误差限eps ，均由用户设定）； 2. 通过数值实验说明，求此线性方程组的近似解时，高斯-赛德尔迭代法的收敛速度比雅可比迭代法的收敛速度要快一些。（用同样精度要求的条件来比较迭代次数） 二、 实验方案(程序源文件) 运用MATLAB 软件编辑M 文件如下： function EX() a=input('请输入系数矩阵a ：'); b=input('请输入矩阵b:'); N=input('请输入最大迭代次数N ：'); esp=input('请输入近似解的误差限：'); if any(diag(a))==0 error('系数矩阵错误，迭代终止！') end D=diag(diag(a)); X0=zeros(size(b)); x1=0; x2=0; x3=0; X1=[x1;x2;x3]; h=inv(D)*b; B=inv(D)*(D-a); B1=triu(B); B2=tril(B); k=1; fprintf('高斯-赛德尔迭代法 \n'); fprintf('第0次迭代得：') disp(X1'); while k<=N x1=h(1,1)+B1(1,:)*X0; X1=[x1;x2;x3]; x2=h(2,1)+B1(2,:)*X0+B2(2,:)*X1; X1=[x1;x2;x3]; x3=h(3,1)+B2(3,:)*X1; X1=[x1;x2;x3]; if norm(X1-X0,inf) 非线性回归预测法——高斯牛顿法(詹学朋) 非线性回归预测法 前面所研究的回归模型，我们假定自变量与因变量之间的关系是线性的，但社会经济现象是极其复杂的，有时各因素之间的关系不一定是线性的，而可能存在某种非线性关系，这时，就必须建立非线性回归模型。 一、非线性回归模型的概念及其分类 非线性回归模型，是指用于经济预测的模型是曲线型的。常见的非线性回归模型有下列几种： (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)泰勒级数展开式。设非线性回归模型为： i=1,2,…,n (3-68) 其中r为待估回归系数，误差项～N(0, )，设： ，为待估回归系数的初始值，将(3-68)式在g点附近作泰勒展开，并略去非线性回归模型的二阶及二阶以上的偏导数项，得 (3-69) 将(3-69)式代入(3-68)式，则 移项： 令： 则：i=1,2，…，n 用矩阵形式表示，上式则为：(3-70) 其中： (3)估计修正因子。用最小平方法对(3-70)式估计修正因子B， 则：(3-71) 设g为第一次迭代值，则： (4)精确度的检验。设残差平方和为： ，S为重复迭代次数，对于给定的允许误差率K，当时，则停止迭代；否则，对(3-71)式作下一次迭代。 (5)重复迭代。重复(3-71)式，当重复迭代S次时，则有：修正因子： 第(S+1)次迭代值： 四、应用举例 设12个同类企业的月产量与单位成本的资料如下表： 表3-9 间接代换法计算表 企业编号单位产品成 本（元） 月产量 1 2 3 4 5 6 7 8 9 10 11 12 160 151 114 128 85 91 75 76 66 60 61 60 10 16 20 25 31 36 40 45 51 56 60 65 (注：资料来源《社会经济统计学原理教科书》第435页) 试配合适当的回归模型分析月产量与单位产品成本之间的关系。 解：(1)回归模型与初始值的选择。根据资料散点图的识别，本数据应配合指数模型：对指数模型两边取对数，化指数模型为线性回归模型，然后施行最小平方法求出初始 值。即： 则上述指数模型变为： 对分别求反对数，得，带入原模型， 得回归模型： 高斯—牛顿迭代法 初始回归模型： 作业：① 分别用J 法和G-S 法求解下列方程，并讨论结果。 123122*********x x x -?????? ??? ?= ??? ? ??? ??????? #include 第十三章 药动学数据的曲线拟合以及常用软件 计算药动学参数，是药代动力学进一步应用的基础。如何测定有关的动力学参数呢？常用的方法是：首先在用药后的若干不同时间，采取血样（或尿样），测定其血药浓度值或尿中药量（这些数值称为实测值或观察值，用C i 表示），这样就有了药物浓度经时曲线数据；然后，依据半对数坐标图，选定一种模型方程 （是时间t 的曲线函数）计算理论估算值(用i C )表示)，按照观察值和理论估算值之 差的平方和（即残差平方和）或加权残差平方和(均用Re 表示)最小的原则，采用适当的算法，求出有关的动力学参数。这种方法，在数学上称为曲线拟合(fitting a curve)。由于所采用的线性药代动力学的模型方程是多指数项之和的函数形式，并且是所含动力学参数的非线性函数，所以这种曲线拟合方法称为非线性最小二乘法。 一、最小二乘法的一般原理 设y 是变量x 的函数，含有m 个待定参数a 1，a 2，…，a m 。记为 y ＝f （x ；a 1，a 2，…，a m ） 若对x 和y 作n 次观察，测得观察值(x 1，y 1)，(x 2，y 2)，…，(x n ，y n )。根据这样一组二维数据，即平面上的若干点，要求确定这个一元函数y ＝f （x ；a 1，a 2，…，a m ）；(i=1，2，…，n)，即一条曲线，使这些点与曲线总体来说尽量接近。并使y 的观察值y i 与理论估算值=i y )f （x i ；a 1，a 2，…，a m ） ；(i=1，2，…，m)的误差平方和，即残差平方和21 ()n e i i i y y ==?∑)2121((,,,,)n i i m i y f x a a a ==?…∑R )i 取得最小值， 或者加权残差平方和21()/n e i i i R y y w ==?∑)2121 ((,,,,))/n i i m i i y f x a a a w ==?…∑ 取得 最小值。其中w i 称为权重系数，在后面的段落会详细讲解。这时Re 有时候也称为目标函数。这就是数据拟合成曲线的思想，简称为曲线拟合。曲线拟合的目的是根据实验获得的数据去建立因变量和自变量之间有效的函数关系，这个函数关系对于药动学来讲就是通过房室模型推导出来的药时曲线公式。根据观察值求出待 311 1.Jacobi迭代法 例1 用jacobi迭代法求解代数线性代数方程组,保留四位有效数字(err=1e-4) 其中A=[8 -1 1;2 10 -1;1 1 -5];b=[1 ;4; 3]。 解:编写jacobi迭代法的函数文件,保存为jacobi.m function [x,k]=jacobi(A,b,x0,eps,N) % 求解Ax=b;x0为初始列向量;eps为误差容限;N为最大迭代次数 % 输出x为近似解;k为迭代次数 n=length(A); x=zeros(n,1); for k=1:N for i=1:n ――――――― end if norm(x-x0,inf) end x0=x; end 编写主程序如下 format long clear A=[8 -1 1;2 10 -1;1 1 -5]; b=[1 ;4; 3]; x0=[0.125; 0.4 ;-0.6 ]; % x0为初始列向量N为最大迭代次数err=1e-4; % err为误差容限 N=25; % N为最大迭代次数 [x,k]=jacobi(A,b,x0,err,N) 得到结果如下 x = 0.22492315625000 0.30561995000000 -0.49388680000000 k = 6 2.Gauss-seidel迭代法 例2 用Gauss-seidel迭代法求解代数线性代数方程组,保留四位有效数字(err=1e-4) 其中A=[8 -1 1;2 10 -1;1 1 -5];b=[1 ;4; 3]。 解:编写Gauss-seidel迭代法的函数文件,保存为gaus.m function [x,k]=gaus(A,b,x0,eps,N) % 求解Ax=b;x0为初始列向量;eps为误差容限;N为最大迭代次数% 输出x为近似解;k为迭代次数 n=length(A); x=zeros(n,1); for k=1:N for i=1:n ―――――― end if norm(x-x0,inf)高斯赛德尔法潮流计算
非线性回归预测法——高斯牛顿法(詹学朋)
高斯-赛德尔迭代法解线性方程组
数值分析练习第五套
数值分析编程及运行结果(高斯顺序消元法)
高斯赛德尔与超松弛迭代法
数值分析版试题及答案
雅可比迭代法与矩阵的特征值
高斯-赛德尔迭代法
非线性回归预测法——高斯牛顿法(詹学朋)知识分享
高斯—牛顿迭代法
雅克比迭代法和高斯-赛德尔法解线性方程组(C++)
第十三章 曲线拟合以及常用的药代动力学软件
matlab Jacobi迭代法Gauss-seidel和SOR迭代