数值分析——线性方程组求解-Toy模板网

这篇具有很好参考价值的文章主要介绍了数值分析——线性方程组求解。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

缘起

清理磁盘的时候偶然发现大二下数值分析的实验作业还在，本着在丢弃之前可以放在网上以备不时之需的原则，我便发了上来。

题目1

分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求解下列线性方程组AX = b，其中A为五对角矩阵(n=20)，b是除第一个分量是1外，其他分量都是0的列向量。
数值分析——线性方程组求解
（1）直接法：观察到该题是形式特殊的五对角矩阵，利用matlab内置的用于产生稀疏矩阵的spdiags函数可以很容易生成系数矩阵A， A的产生代码如下：

e=ones(20,1);
A=spdiags([e -2*e 5*e -2*e e],-2:2,20,20);

利用追赶法(Thomas)可以较合适的解决问题，如果直接利用matlab内置的LU分解函数，显然是毫不费力的。

[L U]=lu(A);

利用上（下）对角系数矩阵在解决线性方程时的优良性质，可以在计算机上轻易得到结果：

b=zeros(20,1);
b(1)=1;
x=U\(L\b);

输出结果为（为了文本美观，输出了x的转置）：

x =

  1 至 11 列
    0.2421    0.0944   -0.0218   -0.0314   -0.0069    0.0049    0.0036    0.0002   -0.0008   -0.0004    0.0001   
  12 至 20 列
    0.0001    0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000

或者为了更好地理解LU分解的细节，我们可以自行编写LU代码，如下：

function [x]=LU(A,b)
 
n=length(b);
U=zeros(n);
L=eye(n);
 
for i=1:n
    U(1,i)=A(1,i);
end
 
for i=1:(n-1)
    for k=(i+1):n
           L(k,i)=(A(k,i)-sum(L(k,1:(i-1))*U(1:(i-1),i)))/U(i,i);
    end
    for k=(i+1):n
        U(i+1,k)=A(i+1,k)-sum(L(i+1,1:i)*U(1:i,k));
    end
end
 
x=zeros(n,1);
y=zeros(n,1);
y(1)=b(1);
 
for i=2:n
    s=0;
    for r=1:(i-1)
        s=s+L(i,r)*y(r);
        y(i)=b(i)-s;
    end
end
 
x(n)=y(n)/U(n,n);
for i=(n-1):-1:1
    s1=0;
    for r=(i+1):n
        s1=s1+U(i,r)*x(r);
        x(i)=(y(i)-s1)/U(i,i);
    end
end

经检验，运行 LU 函数可得相同结果
（2）Jacobi迭代法：

function [ ]=jacobi(A,b,x,iter_max)
 
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=D\(L+U);
f=D\b;
 
y=B*x+f;
r=@(x)norm(A*x-b);
n=1;
x_axis=[n];
y_axis=[r(x)];
z_axis=[norm(x)];
 
while norm(y-x)>=1.0e-6 && n<iter_max
  % fprintf('%f\n',r(x));
    x=y;
    x_axis=[x_axis n+1];
    y_axis=[y_axis r(x)];
    z_axis=[z_axis norm(x)];
    y=B*x+f;
    n=n+1;
end
fprintf('\n共迭代了%d次\n',n);
subplot(1,2,1);
plot(x_axis,y_axis);
xlabel('第i次迭代:(i=1,2,3...)');
ylabel('第i次迭代的残量:(i=1,2,3...)');
subplot(1,2,2);
plot(x_axis,z_axis);
xlabel('第i次迭代:(i=1,2,3...)');
ylabel('第i次迭代x的量:(i=1,2,3...)');
end

在实验中，我选取了零向量作为初始迭代点，A和b的产生与（1）中一模一样，这里不再叙述。iter_max是最大迭代步，当异常发生时，可以及时停止代码而不至于陷入死循环。为了清楚看到迭代的具体速度，这里进行了残量的可视化，输出结果如下：

输入：
>> e=ones(20,1);
>> A=spdiags([e -2*e 5*e -2*e e],-2:2,20,20);
>> b=zeros(20,1);
>> b(1)=1;
>> x=ones(20,1);
>> jacobi(A,b,x,100)
	
输出：
共迭代了100次

数值分析——线性方程组求解
Jacobi法在这个具体问题中不收敛，这是一个意料之中的结果，要使用这种方法，前提是系数矩阵有良好的性质，使得迭代矩阵B的谱半径小于1，我们来算一算该题A的谱半径为多少：

>> D=diag(diag(A));
>> U=-triu(A,1);
>> L=-tril(A,-1);
>> B=D\(L+U);
>> max(abs(eig(B)))

ans =

    1.1743

根据jacobi法收敛的充要条件，我们知道大于1的话迭代将不收敛。顺便地，分享一下，我之前测试用时写地用于产生合适矩阵的函数：

%%
%   BG用于判断谱半径
%   size生成相应大小矩阵
%%
function A=jacobi_matrix(size)
BG=@(x)max(abs(eig(eye(size)-diag(x)\x)));
n=1;
while(n<=500)
    n=n+1;
    A=diag(diag(randn(size)));
    [Q,R]=qr(randn(size));
    a=max(max(abs(A)));
    a=1/(a+0.01);
    A=A*a;
    B=Q*A*Q';
    D=inv(diag(diag(B)));
    L=-triu(B,1);
    U=-tril(B,-1);
    A=D-L-U;
    if(BG(A)<=1.01)
        fprintf('\n success\n');
        return;
    end
end
fprintf('\nfail\n');
end

（3）Gauss_Seidel迭代法：（matlab代码如下：）

%%
%   function name:Gauss_Seidel
%%
function []=Gauss_Seidel(A,b,x,iter_max)
 
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x+f;
r=@(x)norm(A*x-b);
n=1;
x_axis=[n];
y_axis=[r(x)];
z_axis=[norm(x)];
 
while norm(y-x)>=1.0e-6 && n<iter_max
  % fprintf('%f\n',r(x));
    x=y;
    x_axis=[x_axis n+1];
    y_axis=[y_axis r(x)];
    z_axis=[z_axis norm(x)];
    y=B*x+f;
    n=n+1;
end
fprintf('\n共迭代了%d次\n',n);
subplot(1,2,1);
plot(x_axis,y_axis);
xlabel('第i次迭代:(i=1,2,3...)');
ylabel('第i次迭代的残量:(i=1,2,3...)');
subplot(1,2,2);
plot(x_axis,z_axis);
xlabel('第i次迭代:(i=1,2,3...)');
ylabel('第i次迭代x的量:(i=1,2,3...)');
     
end

在运行出结果之前，我感觉大概率也不收敛，结果打脸了，只能说题目出得好啊，这也印证了对于同一个线性方程Ax=b，jacobi法与Gauss-Seidel法的收敛性之间没有必然的联系。

e=ones(20,1);
A=spdiags([e -2*e 5*e -2*e e],-2:2,20,20);
b=zeros(20,1);
b(1)=1;
x=ones(20,1);
Gauss_Seidel(A,b,x,500)

共迭代了19次
数值分析——线性方程组求解

x =

1 至 11 列

    0.2421    0.0944   -0.0218   -0.0314   -0.0069    0.0049    0.0036    0.0002   -0.0008   -0.0004    0.0001

  12 至 20 列

0.0001    0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000

让我们顺便算一下迭代矩阵的谱半径大小，反正又不用我算

>> D=diag(diag(A));
>> U=-triu(A,1);
>> L=-tril(A,-1);
>> B=(D-L)\U;
>> max(abs(eig(B)))

ans =

0.4501

题目2

求 $A x = b$ 的解向量 $x$ ，其中
$A=\begin{pmatrix} 1&2&9\\ 3000&2000&1000\\ 4/10^6 &3/10^6 &2/10^6 \end{pmatrix} ,b=\begin{pmatrix} 1\\ 2000\\ 3/10^6 \end{pmatrix}$
（1）求出系数矩阵的条件数及解向量 $x$
（2）将 $a_{33}$ 改为 $3/10^6$ ， $b_3$ 改为 $4/10^6$ ，求解向量 $\tilde x$ .
（3）令 $P=diag(1,10^{-3},10^6)$ 求解 $P A X = P b$ ，并求系数矩阵 $P A$ 的条件数.
（4）对 $P A$ 中的 $a_{33}$ 和 $P B$ 中的 $b_3$ 给以 $10^{-6}$ 的扰动，求解向量 $\hat x$
结合上述计算结果，讨论以上方程组的性态.
解：
(1)

>> format long
>> A=[1 2 9;3000 2000 1000;4/(1e6) 3/(1e6) 2/(1e6)];
>> cond(A)

ans =

     1.920070697116523e+10

>> condest(A)

ans =

     2.300766669733331e+10

>> b=[1;2000;3/(1e6)];
>> x=A\b

x =

  -0.166666666666667
   1.333333333333333
  -0.166666666666667

(2)

>> A(3,3)=3/(1e6);
>> b(3)=4/(1e6);
>> x=A\b
x =
  -9.499999999999981
  16.499999999999968
  -2.499999999999995

（3）

>> P=diag([1,1e-3,1e6]);
>> x=(P*A)\(P*b)
x =
  -9.499999999999998
  16.499999999999996
  -2.500000000000000
>> cond(P*A)
ans =
     2.655782995863286e+02

>> condest(P*A)
ans =
     3.834999999999999e+02

（4）

>> A(3,3)=4/(1e6);
>> b(3)=5/(1e6);
>> x=(P*A)\(P*b)

x =
  18.500000000000004
 -29.000000000000007
   4.500000000000001

由上述结果可以看到对于病态矩阵（条件数过大），微小的扰动会使解向量的误差变得十分大，其相对误差会变得很大，求出来的解向量没有参考价值。而自行添加一个矩阵P，P起到平衡数量级的作用，可大大降低条件数的大小，虽然求出来的结果还是不精确，但其近似值与精确值已在同个数量级上，这样的解是有意义的，可为接下来的进一步精确提供信息。文章来源地址https://www.toymoban.com/news/detail-445676.html

到了这里，关于数值分析——线性方程组求解的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！