線性方程組AX=B的數(shù)值計(jì)算方法實(shí)驗(yàn)08746

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

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

3、角形線性方程組有如下形式：構(gòu)造一個(gè)程序求解三角形線性方程組?？杉俣ú恍枰儞Q。而且可用第k行消去第k+1行的。核心代碼：#include#include#include#define N 4/矩陣階數(shù)void ColPivot(double cNN+1,double ); /函數(shù)聲明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系數(shù)矩陣為: n;for(i=0;iN;i+)for(j=0;jN;j+)coutsetw(10)cij;coute

4、ndl;cout右側(cè)矩陣 y 為: n;for(i=0;iN;i+)coutsetw(10)ciN;coutendl;cout-endl;ColPivot(c,x); /調(diào)用函數(shù)，進(jìn)行高斯消去法變換cout變換后得到的三角矩陣: n;for(i=0;iN;i+)for(j=0;jN;j+)coutsetw(10)cij;coutendl;cout變換后的右側(cè)矩陣 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; /選出主元后進(jìn)行交換ckj=p;for(j=i+1;jN;j+)p=cji/cii;for(k=i;k=N;k+)cjk-=p*cik; /高斯消去，進(jìn)行計(jì)算for(i=0;i=0;i-) /利用回代法求最終解for(j=N-1;ji;j-)ti-=cij*xj;xi=t

6、i/cii;運(yùn)行結(jié)果：（以具體方程組為例）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(請(qǐng)輸入矩陣元次：n); scanf(%d,&n); printf(請(qǐng)輸入矩陣各項(xiàng)：n); for(i=1;i=n;+i) for(j=1;j=n;+j) scanf(%lf,&aij); printf(請(qǐng)輸入方程組的常數(shù)項(xiàng)：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是否成立，若成立，執(zhí)行步驟（6）；若不成立，執(zhí)行步驟（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是否成立，若成立，執(zhí)行步驟（12）；若不成立，k=p+1，執(zhí)行步驟（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、立，執(zhí)行步驟（12）；若不成立，k=2, Y(1)=E(t(1),J);執(zhí)行步驟（8）判斷kN是否成立，若成立，執(zhí)行步驟（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是否成立，若成立，執(zhí)行步驟（11）；若不成立，X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N)/A(k,k);k=k-1;執(zhí)行步驟（10）.m(1:N,J)=X(1:N);返回步驟（7）.輸出結(jié)果。核心代碼：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