一，算法原理

共轭梯度法可以看作是特殊的迭代法，有迭代法的格式，即首先给出x(0),再由迭代格式

x(k+1)=x(k)+αkd(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)+αkd(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]Td(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)+αAd(k))Td(k)=0{{(-{{r}^{(k)}}+\alpha A{{d}^{(k)}})}^{T}}{{d}^{(k)}}=0(−r(k)+αAd(k))Td(k)=0

其中r(k)=b−Ax(k){{r}^{(k)}}=b-A{{x}^{(k)}}r(k)=b−Ax(k),由此解得最佳步长：

α(k)=r(k)Td(k)d(k)TAd(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−Ax(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)+β0d(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)Ad(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)TAd(k)d(k)TAd(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)+βkd(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−Ax(0)αk=r(k)Td(k)d(k)TAd(k)x(k+1)=x(k)+αkd(k)r(k+1)=b−Ax(k+1)βk=−r(k+1)TAd(k)d(k)TAd(k)d(k+1)=r(k+1)+βkd(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:nalpha=(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)<=e0break;end
end

2，建立课本计算实习3.2方程组的系数矩阵及右端项

%建立课本计算实习3.2方程组的系数矩阵及右端项
function [A,b]=build(n)
A=zeros(n,n);
b=zeros(n,1);
b(1)=-1;
b(n)=-1;
for i=1:nA(i,i)=-2;
end
for j=1:n-1A(j,j+1)=1;
end
for j=2:nA(j,j-1)=1;
end

四，实例分析

计算实例P113：用共轭梯度法求解线性方程组Ax=b，其中

矩阵A的阶数n分别取为100，200，400，指出计算结果是否可靠。
解：输入以下代码

clc;
clear;
n=input('Please input n:');              %输入矩阵A的阶数
e0=input('Please input the accuracy:');
[A,b]=build(n);                        %建立课本例3.2方程组系数矩阵及右端项
[error,x]=gongetidu(A,b,e0);             %使用共轭梯度法进行迭代求解
q=log(error);
plot(q,'linewidth',1.5)
xlabel('迭代次数');
ylabel('log(error)');
title('共轭梯度法迭代误差变化曲线');

例，n=100时，迭代50次后满足精度要求，其误差曲线如下图所示。

Matlab求解线性方程组（一）共轭梯度法相关推荐

MATLAB求解线性方程组的八种方法
MATLAB求解线性方程组的八种方法求解线性方程分为两种方法–直接法和迭代法常见的方法一共有8种直接法 Gauss消去法 Cholesky分解法迭代法 Jacobi迭代法 Gauss-Seid ...
matlab 共轭,求解线性方程组 - 双共轭梯度法
通过为 bicg 提供用来计算 A*x 和 A'*x 的函数句柄(而非系数矩阵 A)来求解线性方程组. 创建一个非对称三对角矩阵.预览该矩阵. A = gallery('wilk',21) + dia ...
[Matlab]求解线性方程组
转自:http://silencethinking.blog.163.com/blog/static/911490562008928105813169/ AX=B或XA=B 在MATLAB中,求解线性 ...
matlab 求解线性方程组之LU分解
线性代数中的一个核心思想就是矩阵分解,既将一个复杂的矩阵分解为更简单的矩阵的乘积.常见的有如下分解: LU分解:A=LU,A是m×n矩阵,L是m×m下三角矩阵,U是m×n阶梯形矩阵 QR分解: 秩分解 ...
MATLAB基础教程（6）——使用matlab求解线性方程组
目录今日任务: 一般方程: 方程组(目前仅讨论方程个数和未知数个数一样的情况): 额外知识咦,咋跑题了左除和右除今日总结: 今日任务: 在数学中经常遇见的一个问题就是方程求解,特别是线性代数中 ...
matlab求解线性方程组
模电题现在看来是不用matlab解方程不可做了orz 绝望,各种绝望,平时不努力到了期末季就焦虑得不行. 左除法就好 x=A/b; 如果有符号变量,用syms声明一下就好. 越来越懒了orz好吧,解线 ...
用matlab求解线性方程组
对方程组Ax=b 1 rank(A)=rank(A,b)=n时,无论m=n还是m>n 有唯一解 m=n时,即方程个数等于未知数个数 A=[1 -1 1 -2;2 0 -1 4;3 2 1 0;- ...
Matlab求解线性方程组的三种方法（wzl）
A=[3 12 1;12 0 2;0 2 3]; b=[2.36;5.26;2.77]; %第一种方法 A\b %当无穷多解时,会使得解中尽量多0 %可算无解方程组%第二种方法 pinv(A)*b % ...
matlab ax=b求x,Matlab求解线性方程组Ax=b的几种常见方法Matlab求解线性方程组Ax=b的几种常见方法...
例如方程组: 法1:左除法 >> A=[3 1 -1;1 2 4;-1 4 5];b=[3.6;2.1;-1.4]; >> x=A\b x = 1.4818 -0.4606 0 ...

Matlab求解线性方程组（一）共轭梯度法

一，算法原理

二，程序框图

三，源代码

四，实例分析

Matlab求解线性方程组（一）共轭梯度法相关推荐

最新文章

热门文章