Matlab求解线性方程组（一）共轭梯度法 |
您所在的位置：网站首页 › maple求共轭 › Matlab求解线性方程组（一）共轭梯度法 |
Matlab求解线性方程组（一）共轭梯度法
|
一，算法原理
共轭梯度法可以看作是特殊的迭代法，有迭代法的格式，即首先给出x(0),再由迭代格式 x ( k + 1 ) = x ( k ) + α k d ( k ) {{x}^{(k+1)}}={{x}^{(k)}}+{{\alpha }_{k}}{{d}^{(k)}} x(k+1)=x(k)+αkd(k)进行迭代，那么关键就是求出两个因素：方向 d ( k ) {{d}^{(k)}} d(k)和步长 α k {{\alpha }_{k}} αk。 首先确定最佳步长 α k {{\alpha }_{k}} αk，假设 x ( k ) {{x}^{(k)}} x(k)和搜索方向 d ( k ) {{d}^{(k)}} d(k)已给定，那么通过一元函数 ϕ ( α ) = f ( x ( k ) + α k d ( k ) ) \phi (\alpha )=f({{x}^{(k)}}+{{\alpha }_{k}}{{d}^{(k)}}) ϕ(α)=f(x(k)+αkd(k))求极小值。令 ϕ ′ ( α ) = 0 \phi '(\alpha )=0 ϕ′(α)=0 ( α \alpha α 为变量，其他为已知数)得到： [ A ( x ( k ) + α d ( k ) ) − b ] T d ( k ) = 0 {{\left[ A\left( {{x}^{(k)}}+\alpha {{d}^{(k)}} \right)-b \right]}^{T}}{{d}^{(k)}}=0 [A(x(k)+αd(k))−b]Td(k)=0 ( − r ( k ) + α A d ( k ) ) T d ( k ) = 0 {{(-{{r}^{(k)}}+\alpha A{{d}^{(k)}})}^{T}}{{d}^{(k)}}=0 (−r(k)+αAd(k))Td(k)=0其中 r ( k ) = b − A x ( k ) {{r}^{(k)}}=b-A{{x}^{(k)}} r(k)=b−Ax(k),由此解得最佳步长： α ( k ) = r ( k ) T d ( k ) d ( k ) T A d ( k ) {{\alpha }_{(k)}}=\frac{{{r}^{{{(k)}^{T}}}}{{d}^{(k)}}}{{{d}^{{{(k)}^{T}}}}A{{d}^{(k)}}} α(k)=d(k)TAd(k)r(k)Td(k)下面确定 d ( k ) {{d}^{(k)}} d(k)，给定x(0)后，由于负梯度是函数下降最快的方向，故第一次迭代取搜索方向: d ( 0 ) = r ( 0 ) = b − A x ( 0 ) {{d}^{(0)}}={{r}^{(0)}}=b-A{{x}^{(0)}} d(0)=r(0)=b−Ax(0)，由上式可算出x(1)。由x(1)出发的搜索方向不再取r(1)而是取 d ( 1 ) = r ( 1 ) + β 0 d ( 0 ) {{d}^{(1)}}={{r}^{(1)}}+{{\beta }_{0}}{{d}^{(0)}} d(1)=r(1)+β0d(0)。 由于 d ( 0 ) {{d}^{(0)}} d(0)和 d ( 1 ) {{d}^{(1)}} d(1)为共轭向量，满足 ( d ( 1 ) A d ( 0 ) ) = 0 ({{d}^{(1)}}A{{d}^{(0)}})=0 (d(1)Ad(0))=0，可求得 β 0 {{\beta }_{0}} β0并给出 β k {{\beta }_{k}} βk的一般表达式： β k = − r ( k + 1 ) T A d ( k ) d ( k ) T A d ( k ) {{\beta }_{k}}=-\frac{{{r}^{{{(k+1)}^{T}}}}A{{d}^{(k)}}}{{{d}^{{{(k)}^{T}}}}A{{d}^{(k)}}} βk=−d(k)TAd(k)r(k+1)TAd(k)再由 d ( k + 1 ) = r ( k + 1 ) + β k d ( k ) {{d}^{(k+1)}}={{r}^{(k+1)}}+{{\beta }_{k}}{{d}^{(k)}} d(k+1)=r(k+1)+βkd(k)算出 d ( k + 1 ) {{d}^{(k\text{+}1)}} d(k+1)。 综上所述得出共轭梯度法的计算公式： { d ( 0 ) = r ( 0 ) = b − A x ( 0 ) α k = r ( k ) T d ( k ) d ( k ) T A d ( k ) x ( k + 1 ) = x ( k ) + α k d ( k ) r ( k + 1 ) = b − A x ( k + 1 ) β k = − r ( k + 1 ) T A d ( k ) d ( k ) T A d ( k ) d ( k + 1 ) = r ( k + 1 ) + β k d ( k ) \left\{ \begin{array}{l} {d^{(0)}} = {r^{(0)}} = b - A{x^{(0)}}\\ {\alpha _k} = \frac{{{r^{{{(k)}^T}}}{d^{(k)}}}}{{{d^{{{(k)}^T}}}A{d^{(k)}}}}\\ {x^{(k + 1)}} = {x^{(k)}} + {\alpha _k}{d^{(k)}}\\ {r^{(k + 1)}} = b - A{x^{(k + 1)}}\\ {\beta _k} = - \frac{{{r^{{{(k + 1)}^T}}}A{d^{(k)}}}}{{{d^{{{(k)}^T}}}A{d^{(k)}}}}\\ {d^{(k + 1)}} = {r^{(k + 1)}} + {\beta _k}{d^{(k)}} \end{array} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧d(0)=r(0)=b−Ax(0)αk=d(k)TAd(k)r(k)Td(k)x(k+1)=x(k)+αkd(k)r(k+1)=b−Ax(k+1)βk=−d(k)TAd(k)r(k+1)TAd(k)d(k+1)=r(k+1)+βkd(k) 二，程序框图
1，主函数（实现共轭梯度求解方程组） %共轭梯度法求线性方程组 % 'A';系数矩阵 % 'b':右端项 % 'e0':求解精度 function [error,x]=gongetidu(A,b,e0) %给定初始向量x0和计算精度e0,error代表误差。 n=length(b); x=zeros(n,1); r=b-A*x; %r为误差 d=r; %d为搜索方向 i=1; for k=1:n alpha=(r'*d)/(d'*A*d); x=x+alpha*d; r1=b-A*x; error(i)=norm(r1); bt=-(r1'*A*d)/(d'*A*d); d=r1+bt*d; r=r1; i=i+1; if norm(r1) |
【本文地址】