举例:
追赶法
用追赶法求解三对角方程组:
(
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
)
(
x
1
x
2
x
3
x
4
x
5
)
=
(
1
0
0
0
0
)
\begin {pmatrix}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{pmatrix}\begin {pmatrix}x_1\\x_2\\x_3\\x_4\\x_5\\\end{pmatrix}=\begin {pmatrix}1\\0\\0\\0\\0\\\end{pmatrix}
⎝⎜⎜⎜⎜⎛2−1000−12−1000−12−1000−12−1000−12⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛x1x2x3x4x5⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛10000⎠⎟⎟⎟⎟⎞
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
求 精确解、条件数、扰动方程
设有线性方程组
A
X
=
b
AX=b
AX=b,其中:
A
=
(
10
7
8
7
7
5
6
5
8
6
10
9
7
5
9
10
)
,
b
=
(
32
23
33
31
)
A=\begin {pmatrix}10&7&8&7\\7&5&6&5\\8&6&10&9\\7&5&9&10\end{pmatrix} , b=\begin {pmatrix}32\\23\\33\\31\end{pmatrix}
A=⎝⎜⎜⎛1078775658610975910⎠⎟⎟⎞,b=⎝⎜⎜⎛32233331⎠⎟⎟⎞
试解答下列问题:
- 给出该方程组的精确解 X X X;
- 计算条件数 C o n d ∞ \rm Cond_\infin Cond∞ ( A ) (A) (A),并由此判定方程组是否病态;
- 若给系数矩阵
A
A
A以扰动;
δ
A
=
0.01
A
\delta A = 0.01A
δA=0.01A,给向量b以扰动;
δ
b
=
(
0.01
,
−
0.01
,
0.01
,
−
0.01
)
T
\delta b = (0.01,-0.01,0.01,-0.01)^T
δb=(0.01,−0.01,0.01,−0.01)T,试求解扰动方程:
( A + δ A ) X ^ = b + δ b , (A+\delta A)\hat X = b+\delta b, (A+δA)X^=b+δb,
并计算其相对误差: e ^ ≜ ∣ ∣ X − X ^ ∣ ∣ ∞ / ∣ ∣ X ∣ ∣ ∞ \hat e \triangleq ||X-\hat X||_\infin/||X||_\infin e^≜∣∣X−X^∣∣∞/∣∣X∣∣∞
主函数
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