1、數值計算方法實驗報告線性方程組AX=B的數值計算方法實驗學號：2013040207012 姓名：梁哲豪實驗描述在自然科學和工程技術中很多問題的解決常常歸結為解線性代數方程組。例如電學中的網絡問題，船體數學放樣中建立三次樣條函數問題，用最小二乘法求實驗數據的曲線擬合問題，解非線性方程組問題，用差分法或者有限元法解常微分方程，偏微分方程邊值問題等都導致求解線性方程組，而且后面幾種情況常常歸結為求解大型線性方程組。線性代數方面的計算方法就是研究求解線性方程組的一些數值解法與研究計算矩陣的特征值及特征向量的數值方法。關于線性方程組的數值解法一般有兩類：直接法：若在計算過程中沒有舍入誤差，經過有限步算術

2、運算，可求得方程組的精確解的方法。迭代法：用某種極限過程去逐步逼近線性方程組精確解的方法。迭代法具有占存儲單元少，程序設計簡單，原始系數矩陣在迭代過程中不變等優點，但存在收斂性及收斂速度等問題。上三角線性方程組的求解：基本算法：高斯消元法：將原方程組化為三角形方陣的方程組： a（k=1,2,n-1; i=k+1,k+2, ,n ;j=k+1,k+2, ,n+1）由回代過程求得原方程組的解：xxLU分解法：將系數矩陣A轉化為A=L*U，L為單位下三角矩陣，U為普通上三角矩陣，然后通過解方程組l*y=b，u*x=y,來求解x。二、實驗內容1、許多科學應用包含的矩陣帶有很多零。在實際情況中很重要的三

3、角形線性方程組有如下形式：構造一個程序求解三角形線性方程組。可假定不需要變換。而且可用第k行消去第k+1行的。核心代碼：#include#include#include#define N 4/矩陣階數void ColPivot(double cNN+1,double ); /函數聲明void main()int i,j;double xN;double cNN+1=1,3,5,7,1,2,-1,3,5,2,0,0,2,5,3,-2,-6,-3,1,4;cout-endl;cout系數矩陣為: n;for(i=0;iN;i+)for(j=0;jN;j+)coutsetw(10)cij;coute

4、ndl;cout右側矩陣 y 為: n;for(i=0;iN;i+)coutsetw(10)ciN;coutendl;cout-endl;ColPivot(c,x); /調用函數，進行高斯消去法變換cout變換后得到的三角矩陣: n;for(i=0;iN;i+)for(j=0;jN;j+)coutsetw(10)cij;coutendl;cout變換后的右側矩陣 y 為: n;for(i=0;iN;i+)coutsetw(10)ciN;coutendl;cout-endl;cout方程的解為: n;for(i=0;iN;i+)cout xi= xiendl;cout-endl;void Col

5、Pivot(double cNN+1,double x)int i,j,k;double p,max;double tN;for(i=0;i=N-2;i+)max=0;k=i;for(j=i+1;jmax)k=j;max=fabs(cji); /選主元if(k!=i)for(j=i;j=N;j+)p=cij;cij=ckj; /選出主元后進行交換ckj=p;for(j=i+1;jN;j+)p=cji/cii;for(k=i;k=N;k+)cjk-=p*cik; /高斯消去，進行計算for(i=0;i=0;i-) /利用回代法求最終解for(j=N-1;ji;j-)ti-=cij*xj;xi=t

6、i/cii;運行結果：（以具體方程組為例）2、（PA=LU：帶選主元的分解法）求解線性方程組AX=B,其中： A= 13核心代碼：#include #include #define L 30 double aLL,bL,lLL,uLL,xL,yL; int main() int n,i,j,k,r; printf(請輸入矩陣元次：n); scanf(%d,&n); printf(請輸入矩陣各項：n); for(i=1;i=n;+i) for(j=1;j=n;+j) scanf(%lf,&aij); printf(請輸入方程組的常數項：n); for(i=1;i=n;+i) scanf(%lf,

7、&bi); for(i=1;i=n;+i) for(j=1;j=n;+j) lij=0; uij=0.0; for(k=1;k=n;+k) for(j=k;j=n;+j) ukj=akj; for(r=1;rk;+r) ukj-=lkr*urj; for(i=k+1;i=n;+i) lik=aik; for(r=1;rk;+r) lik-=lir*urk; lik/= ukk; lkk=1.0; for(i=1;i=n;+i) yi = bi; for(j=1;j0;-i) xi = yi; for(j=i+1;j=n;+j) xi-=uij*xj; xi/= uii; for(i=1;iN-

8、1是否成立，若成立，執行步驟（6）；若不成立，執行步驟（3）。max1,j=max(abs(A(p:N,p);C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=t(p);t(p)=t(j+p-1);t(j+p-1)=d;判斷A(p,p)=0是否成立，若成立，執行步驟（12）；若不成立，k=p+1，執行步驟（5）。判斷kN是否成立，若成立，返回到步驟（2）；若不成立，mult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);k=k+1;返回到步驟（5）。J=1。判斷JN是否成立，若成

9、立，執行步驟（12）；若不成立，k=2, Y(1)=E(t(1),J);執行步驟（8）判斷kN是否成立，若成立，執行步驟（9）；若不成立，Y(k)=E(t(k),J)-A(k,1:k-1)*Y(1:k-1)，k=k+1,返回步驟（8）。X(N)=Y(N)/A(N,N); k=N-1.判斷k 1是否成立，若成立，執行步驟（11）；若不成立，X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N)/A(k,k);k=k-1;執行步驟（10）.m(1:N,J)=X(1:N);返回步驟（7）.輸出結果。核心代碼：function m=newlufact(A,E)%input-A is an N*N

10、 nonsingular matrix% -E is an N*N eyes%output-m is an N*N invmatrix of A %initialize X, Y, the temporary storage matrix C, and the row %permutation information matrix tN,N=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);t=1:N; for p=1:N-1 %find the privot row for column p max1,j=max(abs(A(p:N,p); %in

11、terchange row p and j C=A(p,:); A(p,:)=A(j+p-1,:); A(j+p-1,:)=C; d=t(p); t(p)=t(j+p-1); t(j+p-1)=d; if A(p,p)=0 A is singular. No unique solution break end %calculate multiplier and place in subdiagonal protion of A for k=p+1:N mult=A(k,p)/A(p,p); A(k,p)=mult; A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); endend%solve for Yfo