工程计算——实战：追赶法&扰动分析

举例：

追赶法

用追赶法求解三对角方程组:
$$\left(\begin{array}{ccccc}2& -1& 0& 0& 0\\ -1& 2& -1& 0& 0\\ 0& -1& 2& -1& 0\\ 0& 0& -1& 2& -1\\ 0& 0& 0& -1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\end{array}\right)$$

a=[-1 -1 -1 -1]';
b=[2 2 2 2 2]';
c=a;
d=[1 0 0 0 0]';
n=length(d);

x=chase(a,b,c,d);

function x=chase(a,b,c,d)
n=length(b);
f(1)=c(1)/b(1);
g(1)=d(1)/b(1);
for i=2:n-1
    h(i)=b(i)-f(i-1)*a(i-1);
    f(i)=c(i)/h(i);
    g(i)=(d(i)-g(i-1)*a(i-1))/h(i);
end
g(n)=(d(n)-g(n-1)*a(n-1))/(b(n)-f(n-1)*a(n-1));
x(n)=g(n);
for i=n-1:-1:1
    x(i)=g(i)-f(i)*x(i-1);
end
x
end

这里出现了一些问题 😦
有这样的warning：

数组索引必须为正整数或逻辑值。

出错 chase (line 13)
    x(i)=g(i)-f(i)*x(i-1);

出错 ChaseMain (line 8)
x=chase(a,b,c,d);
 

应该得到的结果：
 x = 0.83333
     0.66667
     0.50000
     0.16667

求 精确解、条件数、扰动方程

设有线性方程组$AX=b$,其中：
$$A=\left(\begin{array}{cccc}10& 7& 8& 7\\ 7& 5& 6& 5\\ 8& 6& 10& 9\\ 7& 5& 9& 10\end{array}\right),b=\left(\begin{array}{c}32\\ 23\\ 33\\ 31\end{array}\right)$$
试解答下列问题：

• 给出该方程组的精确解$X$;
• 计算条件数$\mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{d}}_{\infty }$$(A)$,并由此判定方程组是否病态；
• 若给系数矩阵$A$以扰动；$\delta A=0.01A$,给向量b以扰动；$\delta b=(0.01,-0.01,0.01,-0.01{)}^T$,试求解扰动方程：
  $$(A+\delta A)\hat{X}=b+\delta b,$$
  并计算其相对误差：$\hat{e}\triangleq ||X-\hat{X}|{|}_{\infty }/||X|{|}_{\infty }$
  主函数

A=[10 7 8 7
    7 5 6 5
    8 6 10 9
    7 5 9 10];
b=[32
   23
   33
   31];
%利用Doolittle分解法求该线性方程组精确解
x1=LUX(A,b);
%求解该方程条件数
n=size(A,1);
L=zeros(n);
U=zeros(n);
B=zeros(n);
I=eye(n);

[L,U]=Matrix_LU(A,n);

for i=1:n
    B(:,i)=MatrixInverse(L,U,I(:,i),n);
end
p=input('请输入条件数对应的范数指标（“1”、“2”、“inf”）\n');
mc=MatrixCondition(A,B,p);
fprintf('该线性方程组对应的条件数为：%e',mc);
%加上扰动后利用Doolittle分解法进行求解
A=A+0.01.*A;
delta_b=[0.01
         -0.01
         0.01
         -0.01];
b=b+delta_b;
x2=LUX(A,b);
y1 = mnorminf(x1-x2);
y2 = mnorminf(x1);
e=y1/y2;
e

LUX(A,b):

function x=LUX(A,b)
[n,n]=size(A);
L=zeros(n);U=zeros(n);%用与A相同维度的零矩阵为L、U赋初值
x=zeros(n,1);y=zeros(n,1);
for r=1:n
    for i=r:n
        U(r,i)=A(r,i)-sum(L(r,1:r-1).*U(1:r-1,i)');
        L(i,r)=(A(i,r)-sum(L(i,1:r-1).*U(1:r-1,r)'))/U(r,r);
    end
end

for i=1:n
    y(i)=b(i)-sum(L(i,1:i-1).*y(1:i-1)');
end
L,U
for j=n:-1:1
    x(j)=(y(j)-sum(U(j,j+1:n).*x(j+1:n)'))/U(j,j);
end
x
end

Matrix_LU(A,n):

function [L,U]=Matrix_LU(A,n)

L=zeros(n);U=zeros(n);

for r=1:n
    for i=r:n
        U(r,i)=A(r,i)-sum(L(r,1:r-1).*U(1:r-1,i)');
        L(i,r)=(A(i,r)-sum(L(i,1:r-1).*U(1:r-1,r)'))/U(r,r);
    end
end

end

MatrixInverse(L,U,I(:,i),n)

function x=MatrixInverse(L,U,b,n)

x=zeros(n,1);
y=zeros(n,1);

for i=1:n
    y(i)=b(i)-sum(L(i,1:i-1).*y(1:i-1)');
end

for j=n:-1:1
    x(j)=(y(j)-sum(U(j,j+1:n).*x(j+1:n)'))/U(j,j);
end

end

MatrixCondition(A,B,p)

function mc=MatrixCondition(A,B,p)

switch p
    case 1
        mc = mnorm1(A)*mnorm1(B);
    case 2
        mc = mnorm2(A)*mnorm2(B);
    case inf
        mc = mnorminf(A)*mnorminf(B);
    otherwise
        error('矩阵范数指标错误！')

end

mnorminf:

function y = mnorminf(A)
%矩阵∞范数：最大行绝对和
[m,n] = size(A);
x = zeros(m,1);

for i = 1:m
    for j=1:n
        x(i) = x(i) + abs(A(i,j))
    end
end
y = max(x);
end

最终结果：

%x1
x =

    1.0000
    1.0000
    1.0000
    1.0000
该线性方程组对应的条件数为：4.488000e+03%由于条件数为4488>>1，所以该方程是病态的
L=

    1.0000         0         0         0
    0.7000    1.0000         0         0
    0.8000    4.0000    1.0000         0
    0.7000    1.0000    1.5000    1.0000


U =

   10.1000    7.0700    8.0800    7.0700
         0    0.1010    0.4040    0.1010
         0         0    2.0200    3.0300
         0         0         0    0.5050
%加上扰动后的（即x2）：
x =

    1.8020
   -0.3564
    1.3366
    0.7822

e = 1.3564