用matlab解线性方程组

用matlab解线性方程组

2008-04-12 17:00

一。高斯消去法

1.顺序高斯消去法

直接编写命令文件

a=[]

d=[]'

[n,n]=size(a);

c=n+1

a(:,c)=d;

for k=1:n-1

a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去

end

x=[0 0 0 0]' %回带

x(n)=a(n,c)/a(n,n);

for g=n-1:-1:1

x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)

end

2.列主高斯消去法

*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。

直接编写命令文件

a=[]

d=[] '

[n,n]=size(a);

c=n+1

a(:,c)=d; %(增广)

for k=1:n-1

[r,m]=max(abs(a(k:n,k))); %选主

m=m+k-1; %(修正操作行的值)

if(a(m,k)~=0)

if(m~=k)

a([k m],:)=a([m k],:); %换行

end

a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去end

end

x=[0 0 0 0]' %回带

x(n)=a(n,c)/a(n,n);

for g=n-1:-1:1

x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)

end

3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1]

顺序高斯消去法：提示“Warning: Divide by zero.” x =NaN NaN NaN NaN 列主高斯消去法：x =-1.2460 2.8796 5.5206 -4.3069

由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。

4. 将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“1

5.5”改为

“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。

x =-0.8081 1.8226 3.5568 -2.7047

很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。

这是因为系数矩阵的条件数很大：cond(a)=2.0695e+003

二。迭代法

J迭代公式

a=[5 2 1;-1 4 2;2 -3 10]

d=[-12;20;3]

x=[0;0;0]; %初始向量

stop=1.0e-4 %迭代的精度

L=-tril(a,-1)

U=-triu(a,1)

D=inv(diag(diag(a)))

X=D*(L+U)*x+D*d; % J迭代公式

n=1;

while norm(X-x,inf)>=stop % 时迭代中止否则继续

x=X;

X=D*(L+U)*x+D*d;

n=n+1;

end

X

n

G-S迭代公式

a=[5 2 1;-1 4 2;2 -3 10]

x=[0;0;0]; %初始向量

d=[-12;20;3]

stop=1.0e-4

L=-tril(a,-1)

U=-triu(a,1)

D=(diag(diag(a)))

X=inv(D-L)*U*x+inv(D-L)*d; % G-S迭代公式

n=1;

while norm(X-x,inf)>= stop % 时迭代中止否则继续

x=X;

X=inv(D-L)*U*x+inv(D-L)*d;

n=n+1;

end

X

n

SOR迭代公式

a=[5 2 1;-1 4 2;2 -3 10]

x=[0;0;0]; %初始向量

d=[-12;20;3]

stop=1.0e-4

w=1.44; %松弛因子

L=-tril(a,-1)

U=-triu(a,1)

D=(diag(diag(a)))

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d % SOR迭代公式

n=1;

while norm(X-x,inf)>=stop % 时迭代中止否则继续

x=X;

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d;

n=n+1;

end

X

n

3. a*x=d a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3]

(1)考察用J、G-S及SOR解此方程组的收敛性.

(2)分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求时迭代中止，观察收敛速度。

(1)J迭代矩阵的谱半径：max(abs(eig(D*(L+U))))= 0.506

G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)))= 0.2000 SOR迭代矩阵的谱半径：

max(abs(eig(inv(D-w*L)*((1-w)*D+w*U))))=0.9756

所以用J、G-S及SOR均收敛（且有）。

(2)

J： X =-4.0000 3.0000 2.0000 n =18

G-S： X =-4.0000 3.0000 2.0000 n =8

SOR()： X =-4.0000 3.0000 2.0000 n =482

可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。