实验六 解方程组的直接法

2016-2017学年第1学期数值分析实验
班级姓名实验日期
实验六解方程组的直接法
一、实验目的
熟悉Mathematics中有关矩阵的函数、掌握求解线性方程组的高斯消去法、三角分解法、追赶法的求解方法和计算机实现。

一、实验内容与过程
1 Mathematics中有关矩阵的运算和函数矩阵输入：Shift+Ctri+C
（1）a.b.c or Dot[a, b, c] 矩阵或向量相乘。

In[1]:= {a1,a2,a3}．{b1,b2,b3}
Out[1]= a1 b1+a2b2+a3b3
In[2]:=
1113
236
4326
135
4218
⎛⎫
⎛⎫ ⎪
∙
⎪ ⎪
⎝⎭ ⎪
⎝⎭
//MatrixForm Out[2]=
38231472
33201261
⎛⎫
⎪
⎝⎭
（2）Inverse[A] 求矩阵A的逆矩阵
In[3]:=
123
[203]
111
B Inverse A
⎛⎫
⎛⎫
⎪
⎪
==-
⎪
⎪
⎪
⎪
-
⎝⎭
⎝⎭
//MatrixForm
Out[3]=
3/135/136/13 5/134/133/13 2/131/134/13 -⎛⎫ ⎪
-
⎪ ⎪
-
⎝⎭
（3）Transpose[list] 求矩阵的转置矩阵
In[4]:=Transpose[
123
456
⎛⎫
⎪
⎝⎭
]//MatrixForm Out[4]=
14
25
36
⎛⎫
⎪
⎪
⎪
⎝⎭
（4）Det[m] 求矩阵的行列式
In[5]:=Det[
2346
7852
2312
3541
⎛⎫
⎪
⎪
⎪
⎪
-
⎝⎭
] Out[5]=-112
（5）Tr[list] 求矩阵的迹
In[7]:= mat={{a,b,c},{d,e,f},{g,h,i}})//MatrixForm
Tr[mat]
Out[7]=a+e+i
（6）Eigenvalues[m] 求矩阵的特征值
In[8]:=
21 {,}[]
31
Eigenvalues
λμ
⎛⎫= ⎪
-
⎝⎭
Out[8]=
11 {(1
22
+
（7）LinearSolve [m , b ] 求解线性方程组m .x =b .
In[9]:= 1213m ⎛⎫
= ⎪⎝⎭
；v={5,8};
s=LinearSolve[m,v]
Out[9]= {-1,3}
（8）用三角分解法求解方程组的解
Step1：lub=LUDecomposition [A ] 矩阵A 三角分解A=LU,并将其存放在变量lub 中； Step2：LUBackSubstitution[lub,b], 求解方程Ly=b 和Ux=y 得出方程组的解。

例求解方程组12115108x x ⎛⎫⎛⎫⎛⎫
= ⎪ ⎪ ⎪⎝⎭⎝⎭
⎝⎭
In[10]:= lud=LUDecomposition[1110A ⎛⎫
= ⎪⎝⎭
]
Out[10]= {{{1,1},{1,-1}},{1,2},1}(表示111011101101A ⎛⎫⎛⎫⎛⎫
== ⎪ ⎪⎪-⎝⎭⎝⎭⎝⎭)
In[11]:= LUBackSubstitution [lud,{5,8}] Out[1]= {8,-3}
问题1、用函数LinearSolve 解线性方程组⎪⎩
⎪
⎨⎧=--=++=++2
33322021321321x x x x x x x x
In[9]:= s=LinearSolve[111223130⎛⎫
⎪
⎪ ⎪--⎝⎭
,{0，3，2}]
Out[9]= （写出实验结果）
问题2、用三角分解法解线性方程组 ⎪⎩⎪
⎨⎧=+-=++--=++3
1032202412
25321
321321x x x x x x x x x
In[11]:= lud=LUD ecom position[m=5211422310⎛⎫ ⎪
- ⎪ ⎪-⎝⎭
]
Out[11]= （写出结果）
In[12]:= LUBackSubstitution [lud,{-12,20,3}] Out[12]= （写出结果）
问题3 列主元消去法解线性方程组123412412341342x x -3x -x 13x x 7x 2x 2x 4x 2x 1x x 5x 5+=⎧⎪++=⎪
⎨-++-=-⎪⎪-+=⎩
A={{2,1,-3,-1},{3,1,0,7},{-1,2,4,-2},{1,0,-1,5}};
b={1,2,-1,5};
{m,n}=Dimensions[A];(*求矩阵的行列数*)
Do[temp=Flatten[Abs[Take[A,{k,m},{k}]]];(*从第k列中取出第k到第m个元素*)
kmax=Position[temp,Max[temp]][[1,1]]+k-1;(*找出第k列中最大元素所处的行*)
If[k≠kmax,{A[[k]],A[[kmax]]}={A[[kmax]],A[[k]]};(*矩阵第k行与第kmax行交换*)
{b[[k]],b[[kmax]]}={b[[kmax]],b[[k]]}];(*向量b第k个元素与第kmax个元素交换*)
If[A[[k,k]]≠0,Do[A[[i,k]]=A[[i,k]]/A[[k,k]];
Do[A[[i,j]]=A[[i,j]]-A[[i,k]]*A[[k,j]],
{j,k+1,n}];
b[[i]]=b[[i]]-A[[i,k]]*b[[k]],{i,k+1,m}],
Print["Error"];Break[]],{k,m-1}](*消元*)
b[[m]]=b[[m]]/A[[m,m]];(*求解x
m
*)
Do[b[[i]]=(b[[i]]-Sum[A[[i,j]]*b[[j]],{j,i+1,m}])/A[[i,i]],{i,m-1,1,-1}];(*回代求解*)
b
方程组的解为：
问题4 追赶法解线性方程组
12
123
234
34
3x2x7
-x3x2x11
x3x2x15
x3x9
+=
⎧
⎪++=
⎪
⎨
-++=⎪
⎪-+=
⎩
b={3,3,3,3};
a={-1,-1,-1};
c={2,2,2};
d={7,11,15,9};
n=Length[b];
Moduleodule[{},Do[c[[i-1]]=c[[i-1]]/b[[i-1]];(*计算*)
b[[i]]=b[[i]]-a[[i-1]]*c[[i-1]],{i,2,n}];
d[[1]]=d[[1]]/b[[1]];
Do[d[[i]]=(d[[i]]-a[[i-1]]*d[[i-1]])/b[[i]],{i,2,n}];
Do[d[[i]]=d[[i]]-c[[i]]*d[[i+1]],{i,n-1,1,-1}];d] 方程组的解为：。