beta = zeros(n,1); r(:,1) = b-A*x(:,1); d(:,1) = r(:,1); figure; xlabel('迭代步数'); ylabel('残量'); hold on for k=1:n alpha(k) = r(:,k)'* r(:,k)/(d(:,k)'*A*d(:,k)); x(:,k+1) = x(:,k)+alpha(k)*d(:,k); r(:,k+1) = b-A*x(:,k+1); if norm(r(:,k+1))<epsilon; break; else beta(k) = norm(r(:,k+1))^2/(norm(r(:,k))^2); d(:,k+1) = r(:,k+1)+beta(k)*d(:,k); plot(k,norm(r(:,k+1)),'.k','markersize',15) hold on pause(0.1) end end grid; xlswrite('d:\data_x_result.xls', x(:,k+1), ''\n 所求线性方程组的解为：\n') fprintf('%.8f\n',x(:,k+1))
21
附录 2 生成系数矩阵、右端项以及阶数的 matlab 程序
附录 2 生成系数矩阵、右端项以及阶数的 matlab 程序
clc;clear; n = input('输入系数矩阵的阶数 n： '); A = zeros(n,n); A(1,1:2) = [-2,1]; A(n,n-1:n) = [1,-2]; for i=2:n-1; A(i,i-1:i+1) = [1,-2,1]; end b = zeros(n,1); b(1) = -1; b(n) = -1; csvwrite('d:\data_A.txt',A); csvwrite('d:\data_b.txt',b); csvwrite('d:\data_n.txt',n);
计算方法上机报告
附录 1 共轭梯度法求解线性方程组的 matlab 程序
clear;clc; aa = input('\n 请选择系数矩阵、右端项以及系数矩阵阶数的输 入方式：\n 从文件中输入数据输入 1，\n 从命令窗口输入数据请输 入 2。\n'); if aa==1 A = load('d:\data_A.txt'); b = load('d:\data_b.txt'); n = load('d:\data_n.txt'); end if aa==2 A = input('\n 输入系数矩阵 A(对称正定):\n'); b = input('\n 输入线性方程组的右端项 b:\n'); n = input('\n 输入系数矩阵的阶数 n：\n'); end epsilon = input('\n 输入计算要求的精度 epsilon:\n'); x(:,1) = rand(n,1); alpha = zeros(n,1); %给定初始的向量
22