數值實驗報告
----------------------個人作業，如果有后輩的作業習題一致，可以參考學習，一起交流，請勿直接copy
一、 實驗目的
- 了解並分析LU分解法的優點；
- 追趕法的應用與其與LU分解法的對比；
- 認識迭代法收斂的含義以及迭代法初值和方程組系數矩陣性質對收斂速度的影響。
二、 實驗題目
三、 實驗原理
l LU分解：
- ·如果對A（0）x = b（0）施行第一次消元后化為A（1）x = b(1),則存在L1，使得
L1A（0）=A(1)，L1b（0）= b（1）
一般地，進行k次消元化后為A（k）x = b(k), 則有
LkA（k-1）=A(k)，Lkb（k-1）= b（k）
重復這一過程，最后得到
Ln-1…L2L1A(0) = A(n-1)
Ln-1…L2L1b(0) = b(n-1)
將上三角形矩陣A(n-1)記為U，則 A=LU ，其中
為下三角矩陣。利用高斯消元法實質上產生了一個將A分解為兩個三角形矩陣相乘的因式分解，稱為A的三角形分解或LU分解。
·矩陣分解不一定采用高斯消元法，以下為直接計算的計算公式：
把增廣矩陣A 采用LU 分解格式，即可得到與原方程同解的方
程組。即 Ax=b ⇒ LUx = b ⇒ Ux= L−1b,由U為上三角矩陣，而后回代即可求出原方程的解。
l 追趕法：
求解Ax = b 等價於解兩個二對角線方程組
Ly = b
Ux =y
自上而下解方程組Ly = b 形象地被稱為“追”。
y1 = b1/l11
yi =bi-lii-1yi-1/lii, i = 2, 3, … ,n
自下而上解方程組Ux = y 形象地被稱為“趕”。
xn=yn
xi =yi-uii+1xi+1, i = n-1, … ,2,1
習慣上，上述求解方法稱為“追趕法”。
l 迭代法：
- ·雅克比迭代雅克比迭代法基本思想與迭代法相同是一種逐次逼近的方法。首先給定一個較粗糙的初值，然后采用迭代公式，進行多次迭代，直到滿足所要求的精度為止。
·高斯-賽德爾迭代法基本思想基本與雅克比迭代法相同
四、 實驗內容
- 求矩陣Ａ的ＬＵ分解的主函數如下：
L = new double[n];
U = new double[n];
for (int i = 0; i < s; i++)
{
for (int j = 0; j < s; j++)
{
if (i == j)
L[i*s+j] = 1;
if (i < j)
L[i*s+j] = 0;
if (i > j)
U[i*s+j] = 0;
U[0*s+j] = a[0*s+j];
L[i*s+0] = a[i*s+0] / U[0*s+0];
}
}
for (int k = 1; k < s; k++)
{
for (int j = k; j < s; j++)
{
tmp = 0;
for (int m = 0; m < k; m++)
{
tmp += L[k*s+m] * U[m*s+j];
}
U[k*s+j] = a[k*s+j] - tmp;
}
for (int i = k+1; i < s; i++)
{
tmp = 0;
for (int m = 0; m < k; m++)
{
tmp += L[i*s+m] * U[m*s+k];
}
L[i*s+k] = ( a[i*s+k] - tmp ) / U[k*s+k];
}
}
――――――ＬＵ矩陣求逆矩陣的代碼如下：
//L矩陣求逆
for(int j=0;j<rank;j++)
{
for(int i=j+1;i<rank;i++)
{
double sum=0.0;
for(int k=j;k<i;k++)
{
sum+=Lower[i][k]*ReverseLower[k][j];
}
ReverseLower[i][j]=-sum/Lower[i][i];
}
}
//U矩陣求逆
for(int j=rank-1;j>=1;j--)
{
for(int i=j-1;i>=0;i--)
{
double sum=0.0;
for(int k=i+1;k<=j;k++)
{
sum+=Upper[i][k]*ReverseUpper[k][j];
}
ReverseUpper[i][j]=-sum/Upper[i][i];
}
}
- ＬＵ分解求值Ａｘ＝ｂ的主函數內容如下：
void DirectLU(double a[N][N+1],double x[]) {
int i,r,k,j;
double s[N],t[N];//={-20,8,14,-0.8};
double max;
for(r=0;r<N;r++)
{
max=0;
j=r;
for(i=r;i<N;i++) //求是s[i]的絕對值,選主元
{
s[i]=a[i][r];
for(k=0;k<r;k++)
s[i]-=a[i][k]*a[k][r];
s[i]=s[i]>0?s[i]:-s[i]; //s[i]取絕對值
if(s[i]>max){
j=i;
max=s[i];
}
}
if(j!=r) //選出的主元所在行j若不是r,則對兩行相應元素進行調換
{
for(i=0;i<N+1;i++)
swap(a[r][i],a[j][i]);
}
for(i=r;i<N+1;i++) //記算第r行的元素
for(k=0;k<r;k++){
a[r][i]-=a[r][k]*a[k][i];
}
for(i=r+1;i<N;i++) //記算第r列的元素
{
for(k=0;k<r;k++)
a[i][r]-=a[i][k]*a[k][r];
a[i][r]/=a[r][r];
}
}
for(i=0;i<N;i++)
t[i]=a[i][N];
for(i=N-1;i>=0;i--) //利用回代法求最終解
{
for(r=N-1;r>i;r--)
t[i]-=a[i][r]*x[r];
x[i]=t[i]/a[i][i];
}
}
- 追趕法的具體實現代碼如下：
int flag=1;
cout<<"輸入矩陣A的階數："<<endl;
cin>>M;
b[0]=2;c[0]=1;d[0]=-7;
for(i=1;i<M-1;i++) //輸入各組數據
{
a[i]=1;b[i]=2;c[i]=1;d[i]=-5;
}
a[M-1]=1;b[M-1]=2;d[M-1]=-5;
for(k=0;k<M-1;k++)
{
b[k+1]=b[k+1]-(a[k+1]/b[k])*c[k];
d[k+1]=d[k+1]-(a[k+1]/b[k])*d[k];
}
if(flag)
{
x[M-1]=d[M-1]/b[M-1];
for(i=M-2;i>=0;i--)
{
s=d[i];
s=s-c[i]*x[i+1];
x[i]=s/b[i];
}
cout<<"此方程的解為："<<endl;
for(i=0;i<M;i++)
{
if(i%2==0) cout<<endl;
cout<<"x["<<i<<"]= "<<setprecision(10)<<x[i]<<" ";
//system("pause");
}
cout<<endl;
- Jacobi迭代法的核心算法如下（針對題目4設計的輸入求值改進算法）：
void Jacobi::Iteration()
{
for(i=0;i<n;i++){
sum=0;
for(j=0;j<n;j++)
if(j!=i)
sum+=A[i][j]*x0[j];
if(A[i][i]==0)
{
cout<<"在迭代過程中，系數矩陣主對角線上的數出現為零的情形，無法繼續迭代，程序終止!"<<endl;
exit(0);
}
x[i]=(b[i]-sum)/A[i][i];
}
Judge();
}
- Gauss_Seidel迭代法的核心算法如下（針對題目4設計的輸入求值改進算法）：
void Gauss_Seidel::Iteration()
{
for(i=0;i<n;i++){
sum=sum1=0;
for(j=0;j<i;j++)
sum+=A[i][j]*x0[j];
for(j=i+1;j<n;j++)
sum1+=A[i][j]*x0[j];
if(A[i][i]==0)
{
cout<<"在迭代過程中，系數矩陣主對角線上的數出現為零的情形，無法繼續迭代，程序終止!"<<endl;
exit(0);
}
x[i]=(b[i]-sum-sum1)/A[i][i];
}
Judge();
}
五、 實驗結果
給定矩陣A與向量b
――求A的LU分解；
－－利用A的LU分解求解下列方程：
－Ax = b;
－A2x = b;
－A3x = b
對iii.，若先求LM = A3，再解LMx = b有何缺點？
－――――――――――――――――――――――――――――――――――――――
當n=10時：
――Ax = b的解為
x =（0 .545454，-0.499999，-3.971573，2.5836，-1.444715，8.344654，-5.188799，2.407111，-3.632159，4.545457）－１
――A2x = b 的解為
x =（0.549586，-0.772726，0.249998，8.174827，-4.631004，3.457071，-2.65427，9.026667，-2.272733，4.958682）－１
――A3x = b的解為
x =（0.6883911，-1.172518，0.6363608，-0.1249981，-1.1815，9.874508，-7.902941，1.136398，-4.752079，6.339223）－１
―――――――――――――――――――――――――――――――――――――――
――利用A的LU分解求解A-1，n值自己選定。
當ｎ＝４時：（先求ＬＵ分解，分別求Ｌ和Ｕ的逆矩陣，相乘得Ａ的逆矩陣）
再用現有的LU分解法解此方程組，並對二者進行比較。
由於ｎ＝３００時的解太多，故將其保存在了“追趕法n=300.txt”中，同樣的，n=100時的值也保存在了“追趕法n=100.txt”中
----選取不同初值x0和不同的b，給定迭代誤差用兩種迭代法計算，觀測得到的迭代向量並分析計算結果給出結論
----取定x0及b，將A的主對角線元素成倍放大，其他不變，用簡單迭代法計算多次，比較收斂速度，分析計算結果並給出結論。
六、 實驗結果分析
1.
LU分解矩陣可以較高斯消元法更快速地求解出線性方程組 Ax = b 的解，特別對於像題目1中這樣的特殊下三角形矩陣，使用LU分解可以得到一個主元為1的下三角形矩陣和一個只有主對角線有統一取值的上三角形矩陣；對於求解使其更加方便，大大減少了高斯消元法的運算量；
同時，也因為這個原因，在求解A 的逆矩陣的時候，先行求出 其LU分解的兩個三角形矩陣的逆矩陣，再相乘得到A的逆矩陣，運算量得到了很大的削減，提高了計算效率。
2.
追趕法的基礎是LU分解，所以仍然保持了LU分解的特性，是一種特殊的LU分解，充分利用了系數矩陣的特點，使分解變得更為簡單，得到對三對角線性方程組的快速解法（隨着n的取值的逐漸變大，這兩種解法的差距也明顯地體現了出來，速度和計算量上，LU分解明顯比追趕法要花費大）；
因為追趕法的特性，在求解三對角線性方程組時，利用追趕法要遠優於LU分解法，可以減少運算量，有效提高計算效率。
3.
由兩種不同的迭代法求解題目4中的三對角線性方程組可得，變換不同的初值x0和不同的b,結果的影響基本呈線性變換，受兩者影響；固定其中某個值不變只改變另一組時，其可以明顯察覺出其中的變化。但是兩者迭代方法收斂的速度相差不大，高斯-賽德爾略優於雅克比迭代法；
當取定x0和b，將Ａ的主對角線元素成倍放大，其他不變時，可以明顯地看出，A的主對角元素越大，迭代法迭代的次數就越少，也就意味着收斂越快；因為主元在這個方程中起到的影響因素明顯變大，其他數值的影響被削弱，更容易找到收斂的結果，迭代次數更少，速度更快；
在題目5中，由②中的特殊結果表明，對於簡單的二元方程組，迭代法的實用性明顯變差，而且特別依賴於初值的設定，如果初值設定的不正確或者稍有差別，結果就會非常復雜，甚至於偏差只有1，卻導致迭代了近200次，結果卻越來越差的計算值。所以，對於簡單的線性方程組，用簡單的高斯消元法反而會更加簡單。
--------------------------------------
因為直接從word實驗報告中copy過來的，排版會有些奇怪，不過一些重點的部分已經用截圖替換了，應該不會影響閱讀
有關的算法程序關鍵部分已經在實驗報告中展示，詳細的代碼，可以自行補充或者評論私信。