实验2 线性方程组的直接解法

实验2 线性方程组解的直接解法

【实验目的】

1、验证高斯消去法和三角分解法

2、掌握直接求解线性方程组的常用算法：列主元高斯消去法、LU 分解法等

3、记录运行结果，回答问题，完成实验报告

【实验要求】

根据题目要求，用Matlab 完成下列实验内容。书写实验报告。

【实验内容】

1、利用Matlab 求解线性方程组

考虑下面给出的线性方程组形式为Ax =B 。其中，A 和B 均为给定矩阵

1111n1,n nn n a a b A B a a b ???? ? ?== ? ? ? ?????

可以由给定的A 和B 矩阵构造出解的判定矩阵C （可以用C=[A B]命令来生成C ）

11111

n n nn n a a b C a a b ?? ?= ? ???

这里先可以给出线性方程组有解的判定定理：

A=[1 2 3 ;2 -1 -1;1 -2 -2;1 -1 -1];B=[2;1;-1;0];

C=[A,b];

C,

rank(A),

rank(C),

x=A\B,

x=[1;2;-1];

norm(A*x-B)

答案是：

C =

1 2 3 2

2 -1 -1 1

1 -

2 -2 -1

1 -1 -1 0

ans =

3

ans =

3

x =

1.0000

2.0000

-1.0000

ans =

（2）当rank(A )=rank(C )=r ＜n 时，方程组有无穷多解。具体求解方法留给同学

们实验课后自学。

（3）当rank(A )＜rank(C )时，方程组为矛盾方程，只能解出方程的最小二乘解。

这部分内容将在第四章学习。

【问题1】：求解下面给出的线性方程组。记录所用命令和求解结果。

123454

3214132434

1322x ???? ? ? ? ?= ? ? ? ?????

输入命令：

A=[1 2 3 4;4 3 2 1;1 3 2 4;4 1 3 2];

B=[5;4;3;2];

C=[A,B];

C,

rank(A),

rank(C),

x=A\B,

x=[-1.8;1.8667;3.8667;-2.1333];

B=[5;4;3;2];

norm(A*x-B)

答案是：

C =

1 2 3 4 5

4 3 2 1 4

1 3

2 4 3

4 1 3 2 2

ans =

4

ans =

4

x =

-1.8000

1.8667

3.8667

-2.1333

ans =

5.0990e-004

【思考1】：Matlab 对线性方程组求解是否要考虑主元的选取？为什么？可以用

下列例子试试

121

20.000112x x x x +=??+=?

答：

A=[0.0001 1;1 1];

B=[1;2];

C=[A,B];

C,

rank(A),

rank(C),

x=A\B,

X=[1.0001;0.9999];

norm(double(A*X-B))

输出答案是：

C =

0.0001 1.0000 1.0000

1.0000 1.0000

2.0000

ans =

2

ans =

2

x =

1.0001

0.9999

ans =

1.0000e-008

2、用Matlab 编写高斯消去法具体实现的函数

要求输入参数包含方程组的系数矩阵及常数项向量，输出方程组的解，并使用问题1的例子进行调试。

高斯顺序消去法:设线性方程组AX ＝b 。

1.消元过程

设0)(≠k kk a ,对1,,2,1-=n k 计算

()()(1)()()(1)()()/,1,2,,k k ik ik kk k k k ij ij ik kj k k k i i ik k m a a a a m a i j k k n b b m b ++?=?=-=++??=-?K 其中

⒉回代过程

()()()()()1/1,,2,1()/n n n n nn n i i i i i ij j ii j i x b a i n x b a x a =+?=?=-?=-??∑K 其中。

参考实验框架

输入：

function [x] = guass(A,b)

%开始计算，先赋初值

A=input('输入系数矩阵A的值');

b=input('输入值矩阵b的值');

n=length(A);

x = zeros(n,1);

%消元过程

for k = 1:n-1

for i=j=k+1:n

m(i:i,k:k)=A(i:i,k:k)/A(k:k,k:k);

A(i:i,j:j)=A(i:i,j:j)-m(i:i,k:k)*A(k:k,j:j);

b(i:i)=b(i:i)-m(i:i,k:k)*b(k:k);

end

end

%回代过程

for k = n:-1:1

for i=j=k+1:n

x(n:n)=b(n:n)/A(n:n,n:n);

x(i:i)=(b(i:i)-sum(A(i:i,i+1:n)*x(i+1:n))/A(i:i,i:i));

end

end

想试下用实例做的，不会做啊。

A=[1 2 3 4;4 3 2 1;1 3 2 4;4 1 3 2];

B=[5;4;3;2];

n=length(A);

x = zeros(n,1);

%消元过程

for k = 1:n-1

i=k+1;

j=j+1;

for i=k+1:n

m=A(i:i,k:k)/A(k:k,k:k);

A(i:i,k:n)=A(i:i,k:n)-m*A(k:k,k:n);

end

end

C=A

%回代过程

for k = n-1:1

i=k+1

x(n:n,n:n)=b(n:n,n:n)/A(n:n,n:n)

x(i:i)=(b(i:i)-sum(A(i:i,i+1:n)*x(i+1:n))/A(i:i,i:i)); end