共轭梯度法、最速下降法求解大规模稀疏方程组【Matlab】

77 阅读 0 评论 51 点赞

我是靠谱客的博主辛勤缘分，最近开发中收集的这篇文章主要介绍共轭梯度法、最速下降法求解大规模稀疏方程组【Matlab】，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

针对此题，可分别用共轭梯度法、最速下降法求解线性方程组。

程序如下：

附录1 共辄梯度法求解大规模稀疏方程组程序

附录2 三对角矩阵A、右端项b生成程序

附录3 最速下降法求解线性方程组程序

% 附录1 共轭梯度法求解大规模稀疏方程组程序
%% 利用共轭梯度法求解大规模稀疏方程组
clear       %清除变量
clc         %清除命令行窗口代码
aa=input('n请选择系数矩阵A、右端项b的输入方式：n从文件中输入数据请输入0，n从命令行窗口输入数据请输入1n');
if aa==0
    A = load('data_A.txt');
    b = load('data_b.txt');
end
if aa==1
    A=input('n请输入系数矩阵A(对称正定)：n');
    b=input('n请输入线性方程组的右端项b：n');
end
[n1,n2]=size(A);           %求解矩阵维数
n3=length(b);
if n1==n3 && n2==n3        %验证A、b维数是否相等
    n=n1;
else
    fprintf('nA、b维数不相等，程序结束。n')
    return;
end
[R,flag] = chol(A);        %验证A是否对称正定
if flag ~= 0
    fprintf('nA不是对称正定矩阵，程序结束。n');
    return;
end
accuracy = input('n请输入计算要求的精度：n');   %给定计算要求的精度
x(:,1) = zeros(n,1);                               %给定初始向量
r(:,1) = b-A*x(:,1);
y1(1)=log10(norm(r(:,1)));                         %初始残差取对数
d(:,1)= r(:,1);                                     %搜索方向
for k=1:n
    alpha(k)=(r(:,k))'*d(:,k)/((d(:,k))'*A*d(:,k));  %最佳步长
    x(:,k+1)= x(:,k)+alpha(k)*d(:,k);
    r(:,k+1)= b-A*x(:,k+1);
    y1(k+1)=log10(norm(r(:,k+1)));       %纵轴为误差取对数
    if norm(r(:,k+1))<=accuracy
        break;
    else
        beta(k)=-(r(:,k+1))'*A*d(:,k)/((d(:,k))'*A*d(:,k));
        d(:,k+1)=r(:,k+1)+beta(k)*d(:,k);
    end
end
figure;
plot([0:k],y1,'- r o', 'markersize',2)
xlabel('迭代步数');
ylabel('误差(取对数)');
title(['n=',num2str(n),'时的收敛速度图']);
grid;
xlswrite('data_x_result.xls', x(:,k+1), 'sheet1');   
%将x结果以“.xls”格式保存至“程序所在文件夹”目录下
fprintf('n所求线性方程组的解为：n')          %显示x数值解结果
fprintf('%.6fn',x(:,k+1))

% 附录2 三对角矩阵A、右端项b生成程序
%% 三对角矩阵A、右端项b生成
clear       %清除变量
clc         %清除命令行窗口代码
n = input('请输入系数矩阵A的阶数n：n');
a = input('请输入系数矩阵A第n行第n-1、n、n+1个数：n');
A= zeros(n,n);
A(1,1:2)=-[a(2),a(3)];
A(n,n-1:n)=-[a(1),a(2)];
for i=2:n-1
    A(i,i-1:i+1)=-[a(1),a(2),a(3)];
end
b = zeros(n,1);
b(1)=1;
b(n)=1;
csvwrite('data_A.txt', A);
csvwrite('data_b.txt', b);

% 附录3 最速下降法求解线性方程组程序
%% 最速下降法，求解线性方程组Ax=b
clear       %清除变量
clc         %清除命令行窗口代码
aa=input('n请选择系数矩阵A、右端项b的输入方式：n从文件中输入数据请输入0，n从命令行窗口输入数据请输入1n');
if aa==0
    A = load('data_A.txt');
    b = load('data_b.txt');
end
if aa==1
    A=input('n请输入系数矩阵A(对称正定)：n');
    b=input('n请输入线性方程组的右端项b：n');
end
[n1,n2]=size(A);           %求解矩阵维数
n3=length(b);
if n1==n3 && n2==n3        %验证A、b维数是否相等
    n=n1;
else
    fprintf('nA、b维数不相等，程序结束。n')
    return;
end
[R,flag] = chol(A);        %验证A是否对称正定
if flag ~= 0
    fprintf('nA不是对称正定矩阵，程序结束。n');
    return;
end
accuracy = input('n请输入计算要求的精度：n');   %给定计算要求的精度
x(:,1) = rand(n,1);                     %给定初始向量
r(:,1)=b-A*x(:,1);   y(1)=log10(norm(r(:,1)));   %初始残差取对数
k=1;
while norm(r(:,k))> accuracy
    alpha=(r(:,k))'*r(:,k)/(r(:,k)'*A*r(:,k));
    x(:,k+1) =x(:,k)+alpha*r(:,k);
    r(:,k+1)=b-A*x(:,k+1);
    y(k+1)=log10(norm(r(:,k+1)));           %残差取对数
    k=k+1;
end
figure;
plot([0:k-1],y,'- r o', 'markersize',2)
xlabel('迭代步数');
ylabel('误差');
title(['n=',num2str(n),'时的收敛速度图']);
grid;
xlswrite('data_x_result.xls', x(:,k), 'sheet1');
fprintf('n所求线性方程组的解为：n')
fprintf('%.6fn',x(:,k))

此处仅展示了代码实现，具体算法原理可参考《数值分析》等有关书籍，或在博主“资源”下载页面，下载相关文档进行查看。

计算方法-上机作业-示例【仅供交流参考】-统计分析文档类资源-CSDN文库https://download.csdn.net/download/m0_56425991/87264676