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

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

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

程序如下：

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

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

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

% 附录1 共轭梯度法求解大规模稀疏方程组程序
%% 利用共轭梯度法求解大规模稀疏方程组
clear       %清除变量
clc         %清除命令行窗口代码
aa=input('n请选择系数矩阵A、右端项b的输入方式：n从文件中输入数据请输入0，n从命令行窗口输入数据请输入1n');
if aa==0A = load('');b = load('');
end
if aa==1A=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;
elsefprintf('nA、b维数不相等，程序结束。n')return;
end
[R,flag] = chol(A);        %验证A是否对称正定
if flag ~= 0fprintf('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:nalpha(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))<=accuracybreak;elsebeta(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-1A(i,i-1:i+1)=-[a(1),a(2),a(3)];
end
b = zeros(n,1);
b(1)=1;
b(n)=1;
csvwrite('', A);
csvwrite('', b);% 附录3 最速下降法求解线性方程组程序
%% 最速下降法，求解线性方程组Ax=b
clear       %清除变量
clc         %清除命令行窗口代码
aa=input('n请选择系数矩阵A、右端项b的输入方式：n从文件中输入数据请输入0，n从命令行窗口输入数据请输入1n');
if aa==0A = load('');b = load('');
end
if aa==1A=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;
elsefprintf('nA、b维数不相等，程序结束。n')return;
end
[R,flag] = chol(A);        %验证A是否对称正定
if flag ~= 0fprintf('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))> accuracyalpha=(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文库