C++_數值分析列主元高斯消去法_列主元LU分解法_Jacobi迭

1、數值分析上機實驗報告華南理工大學 08級計算機科學與技術2班1.三次樣條插值（初值條件1）：P52.9、給定函數的函數表和邊界條件，求三次樣條插值函數，并求的近似值。源代碼：yangtiao.cpp#include#includevoid main()int choice = 0;int n = 2;double xx,*x, *y,*a,*b,*a1,*b1,*h,*m;cout請輸入插值節(jié)點個數n：n;x = new doublen; y = new doublen;a = new doublen; b = new doublen;a1 = new doublen; b1 = new do

2、ublen;h = new doublen-1; m = new doublen+1;cout請輸入n個插值的節(jié)點(xi，yi)：endl;for (int i = 0 ; i xiyi;for (int j = 0 ; j n-1 ; j+)hj = xj+1 - xj;cout請輸入待估點xx：xx;cout請選擇邊界條件：choice;switch(choice)case 1:double temp1,temp2;a0 = 0;an-1 = 1;cout請輸入邊界條件的兩個一階微商值s(x1)與s(xn)：temp1temp2;b0 = 2*temp1;bn-1 = 2*temp2;br

3、eak;case 2:a0 = 1;an-1 = 0;b0 = 3/h0*(y1-y0);bn-1 = 3/hn-2*(yn-1-yn-2);break;for (int k = 1 ; kn-1 ; k+)ak = hk-1/(hk-1+hk);bk = 3*( (1-ak)/hk-1*(yk-yk-1) + ak/hk*(yk+1-yk) ); a10 = -a0/2;b10 = b0/2;for (int l = 1; l=0; j-)mj = a1j * mj+1 + b1j;/判別xx所在區(qū)間并輸出結果coutn插值結果為：;for(k = 0 ;k n-1 ;k+)if (xk x

4、x )double output = 0;output = (1+2*(xx-xk)/(xk+1-xk)*pow(xx-xk+1)/(xk-xk+1),2) *yk+(1+2*(xx-xk+1)/(xk-xk+1)* pow(xx-xk)/(xk+1-xk),2) *yk+1+(xx-xk)* pow(xx-xk+1)/(xk-xk+1),2) *mk+(xx-xk+1)* pow(xx-xk)/(xk+1-xk),2) *mk+1;cout output endl;break;delete x;delete y; delete a; delete b; delete a1;delete b1

5、; delete h; delete m;運行結果截圖：2.三次樣條插值（初值條件2）：P52.10、給定函數的函數表和邊界條件，求三次樣條插值函數，并求的近似值。源代碼：yangtiao.cpp(同上)運行結果截圖：3.自動選取步長梯形法：P977、使用自動選取步長梯形法計算積分的近似值。（給定=0.01）源代碼：SelfSelLength.cpp#include#includedouble fun(double a)return 2/( 1+a*a );double SelfSelLength(double R_a,double R_b,double e)double h = (R_b-R

6、_a)/2;double R1 = (fun(R_a)+fun(R_b) * h;int n = 1;double R0;double S;double E;do /每當誤差值不符合要求時，計算下一個result值R0 = R1;S = 0;for (int k =1 ;k 3*e);return R1;void main()double a,b,e;cout 請依次輸入待求積分函數的下界a、上界b 及 精度要求e: a b e;cout 自動選取步長梯形法可求得積分: SelfSelLength (a,b,e)endl; 運行結果截圖：4.Romberg求積法：P978、使用Romberg求

7、積法計算積分的近似值。（給定=0.01，且?。┰创a：Romberg.cpp#include#includestatic double Tri128128;double fun(double a)return sqrt(a);double Romberg(double R_a,double R_b,double e)Tri00 = (R_b-R_a)/2*(fun(R_a)+fun(R_b);int k = 0;double E;do /每當誤差值不符合要求時，計算下一行的Tri值k+;double temp = 0;/計算T0k的數值for (int i = 1 ; i=pow(2,k-1)

8、 ;i+)temp = temp + fun(R_a + (2*i-1)* (R_b - R_a)/pow(2,k);Tri0k = 0.5 * (Tri0k-1 + (R_b-R_a)/pow(2,k-1) * temp);for (int m = 1 ; m e);return Trik0;void main()double a,b,e;cout 請依次輸入待求積分函數的下界a、上界b 及 精度要求e: a b e;cout 按Romberg求積法可求得積分: Romberg (a,b,e)endl;運行結果截圖：5.列主元高斯消去法：測試矩陣為：=源代碼：Guess_Eliminatio

9、n.cpp/列主元高斯消去法#include#include#define N 4 /矩陣的維數，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系數矩陣 static double BN=1,0,1,0; /右端項static double XN;int i,j,k; /計數器void main()for(k = 0; k N-1 ;k+)/選取最大主元int index = k;for(i = k; i N ;i+)if(fabs(Aindexk) fabs(Aik)index = i;/交換行double t

10、emp;for( i = k ; i N ;i+ )temp = Aindexi;Aindexi = Aki;Aki = temp;temp = Bindex;Bindex = Bk;Bk = temp;for(i = k+1; iN; i+)double T = Aik/Akk;Bi = Bi - T * Bk;for ( j = k+1 ; j =0 ; i-)double Temp = 0;for (int j = i+1; jN ;j+)Temp = Temp + Aij * Xj;Xi = (Bi - Temp) /Aii;cout 線性方程組的解（X1，X2，X3.Xn）為：end

11、l; for( i = 0; i N ;i+)cout Xi ;運行結果截圖：6.列主元LU分解法：測試矩陣為：=源代碼：LU_Decomposition.cpp#include#include#define N 4 /矩陣維數，可自定義static double ANN; /系數矩陣static double BN; /右端項static double YN; /中間項static double XN; /輸出static double SN; /選取列主元的比較器int i,j,k; /計數器void main()cout 請輸入線性方程組（ai1，ai2，ai3.ain, yi）：end

12、l; for ( i = 0; i N ;i+)for (int j = 0; j Aij;cin Bi;for (k = 0 ; k N; k+)/選列主元int index = k;for (i = k ; i N; i+) double temp = 0; for (int m = 0 ; m k ;m+) temp = temp + Aim*Amk; Si = Aik - temp; if(Sindex Si) index = i; /交換行 double temp; for( i = k ; i N ;i+ ) temp = Aindexi; Aindexi = Aki; Aki =

13、 temp; temp = Bindex; Bindex = Bk; Bk = temp; / 構造L、U矩陣 for (j = k ; j N; j+) double temp = 0; for (int m = 0 ;m k; m+ ) temp = temp + Akm * Amj; Akj = Akj - temp; /先構造U一行的向量 for( i = k+1; i N; i+) double temp = 0; for (int m =0 ; m k; m+ ) temp = temp + Aim * Amk; Aik = (Aik - temp)/Akk; /再構造L一列的向量

14、 /求解LY = B Y0 = B0; for (i = 1; i N ; i+) double temp = 0; for (int j =0 ; j = 0 ; i- )double temp = 0;for (int j =i+1 ; j N; j+ )temp = temp + Aij * Xj;Xi = (Yi - temp)/Aii;/打印Xcout 線性方程組的解（X1，X2，X3.Xn）為：endl; for( i = 0; i N ;i+)cout Xi ;運行結果截圖：7.簡單迭代法（Jacobi迭代）：測試矩陣為：= 取精度為esp = 0.001，最大迭代次數M = 1

15、00源代碼：Jacobi.cpp#include#include#define N 4 /矩陣的維數，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系數矩陣 static double BN=1,0,1,0; /右端項static double YN; /輸出比較項static double YN; static double XN; /輸出項static double GN; /X = BX + G的G矩陣int i,j,k; /計數器double eps = 0.001;int M = 100;bool di

16、stance() /求兩輸出項的差的范數是否滿足精度要求double temp = 0;for (i = 0 ;i eps)return false;elsereturn true; /滿足精度要求則結束程序void main()/形成迭代矩陣B，存放到A中for (i = 0 ;i N;i+)if (fabs(Aii) eps)cout 打印失敗endl;return;double T = Aii;for ( j = 0 ; j N;j+)Aij = -Aij/T;Aii = 0;Gi = Bi/T;int counter = 0;while (counter M)/迭代for (i = 0

17、;i N; i+)double temp = 0;for (j = 0;jN; j+)temp = temp + Aij*Yj;Xi = Gi + temp;if (distance()=true)break;else/交換X，Y向量；for( i = 0; i N ;i+)Yi = Xi;counter+;/打印Xcout 迭代次數為：counter次。該線性方程組的解（X1，X2，X3.Xn）為：endl;for( i = 0; i N ;i+)cout Xi ;運行結果截圖：8.Seidel迭代：測試矩陣為：= 取精度為esp = 0.001，最大迭代次數M = 100源代碼：Seide

18、l.cpp#include#include#define N 4 /矩陣的維數，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系數矩陣 static double BN=1,0,1,0; /右端項static double YN; /輸出比較項static double XN; /輸出項static double GN; /X = BX + G的G矩陣int i,j,k; /計數器double eps = 0.001;int M = 100;bool distance() /求兩輸出項的差的范數是否滿足精度要求d

19、ouble temp = 0;for (i = 0 ;i eps)return false;elsereturn true; /滿足精度要求則結束程序void main()/形成迭代矩陣B，存放到A中for (i = 0 ;i N;i+)if (fabs(Aii) eps)cout 打印失敗endl;return;double T = Aii;for ( j = 0 ; j N;j+)Aij = -Aij/T;Aii = 0;Gi = Bi/T;int counter = 0;while (counter M)/迭代for (i = 0;i N; i+)double temp = 0;for

20、(j = 0;jN; j+)temp = temp + Aij*Xj;Xi = Gi + temp;if (distance()=true)break;else/交換X，Y向量；for( i = 0; i N ;i+)Yi = Xi;counter+;cout 迭代次數為：counter次。該線性方程組的解（X1，X2，X3.Xn）為：endl;for( i = 0; i N ;i+)cout Xi ;9.松弛法（SOR迭代）測試矩陣為： 取松弛系數w = 1.46= 取精度為esp = 0.001，最大迭代次數M = 100源代碼：SOR.cpp#include#include#define N 4 /矩陣的維數，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系數矩陣 static double BN = 1,0,1,0; /右端項static double YN; /輸出比較項static double X