典型數(shù)值算法的說明C++語言程序的設(shè)計開發(fā)

1、倏4獵菽人牽數(shù)值方法課程設(shè)計說明書典型數(shù)值算法的C+語言程序設(shè)計學(xué)生姓名：晏瑞學(xué) 號：200912010127院（系）：理學(xué)院專 業(yè)：數(shù)學(xué)與應(yīng)用數(shù)學(xué) 091班指導(dǎo)教師：劉海峰2011 年6月15日陜西科技大學(xué)數(shù)值計算課程設(shè)計任務(wù)書理學(xué)院應(yīng)用數(shù)學(xué)專業(yè)數(shù)學(xué)091班級學(xué)生：晏瑞 題目：典型數(shù)值算法的C+語言程序設(shè)計課程設(shè)計從 2011年5月20日起到 2011年6月25日1、課程設(shè)計的內(nèi)容和要求（包括原始數(shù)據(jù)、技術(shù)要求、工作要求等）：每人需作10個算法的程序、必做6題、自選4題。對每個算法要求用C+語言進行編程。必選題：1、經(jīng)典四階龍格庫塔法解一階微分方程組2、高斯列主元法解線性方程組3、牛頓法解非

2、線性方程組4、龍貝格求積分算法5、三次樣條插值算法（壓緊樣條）用 C+語言進行編程計算依據(jù)計算2S果，用Matlab畫圖并觀察三次樣條插值效果。6、M次多項式曲線擬合，據(jù)計算結(jié)果，用 Matlab畫圖并觀察擬合效果。自選題：自選4道其他數(shù)值算法題目.每道題目重選次數(shù)不得超過5次.2、對課程設(shè)計成果的要求包括圖表、實物等硬件要求：1）提交課程設(shè)計報告按照算法要求，用 C+語言設(shè)計和開發(fā)應(yīng)用程序，提交由算法說明；程序 設(shè)計說明；系統(tǒng)技術(shù)文檔（包括系統(tǒng)各模塊主要流程圖，軟件測試方案與測試記 錄、軟件調(diào)試和修改記錄、測試結(jié)論、運行情況記錄），系統(tǒng)使用說明書，源程序代碼為附錄構(gòu)成的課程設(shè)計報告。陜西科技

3、大學(xué)2）課程設(shè)計報告版式要求打印版面要求：A4紙，頁邊距：上 2cm,下2cm,左2.5cm、右2cm；字 體：正文宋體、小四號；行距：固定值 20;頁眉1.5cm ,頁腳1.75cm；頁碼位 于頁腳居中打?。黄鏀?shù)頁頁眉“數(shù)值計算課程設(shè)計”，偶數(shù)頁頁眉“算法名稱”, 頁眉宋體小5號；段落及層次要求：每節(jié)標(biāo)題以四號黑體左起打?。ǘ吻岸魏蟾?0.5行），節(jié)下為小節(jié)，以小四號黑體左起打?。ǘ吻岸魏蟾?.5行）。換行后以小四號宋體打印正文。節(jié)、小節(jié)分別以1、1.1、1.1.1依次標(biāo)出，空一字符后接各部分的標(biāo)題。當(dāng)論文結(jié)構(gòu)復(fù)雜，小節(jié)以下的標(biāo)題，左起頂格書寫，編號依次用（1）、（2） 或1）、2）順序表示

4、。字體為小四號宋體。對條文內(nèi)容采用分行并敘時，其編號用（a）、（b）或a）、b）順序 表示，如果編號及其后內(nèi)容新起一個段落，則編號前空兩個中文字符。3）設(shè)計報告裝訂順序與規(guī)范封面數(shù)值計算課程設(shè)計任務(wù)書目錄數(shù)值計算設(shè)計課程設(shè)計報告正文設(shè)計體會及今后的改進意見參考文獻（資料）左邊緣裝訂3、課程設(shè)計工作進度計劃：時間設(shè)計任務(wù)及要求第16周編寫和調(diào)試程序并按要求撰寫設(shè)計報告指導(dǎo)教師：日期：教研室主任：日期：目錄1 .經(jīng)典四階龍格庫塔法解一階微分方程 -7 -1.1 算法說明： -71.2 ,算法程序： -81.3 ,運行結(jié)果 -10-2 .高斯列主元法解線性方程組 -10-2.1 算法說明： -11-

5、2.2 算法程序： -11 -2.3 運行結(jié)果： -14-3 .牛頓法解非線性方程組 -16-3.1 算法說明 -16-3.2 算法程序： -16-3.3 運行結(jié)果： -20-4 .龍貝格求積分算法 -22 -4.1 算法說明 -22-4.2 算法程序： -22-4.3 運彳丁結(jié)果： -24 -5 .三次緊壓樣條插值 -24 -5.1 算法說明： -24-5-數(shù)值計算課程設(shè)計5.2 算法程序： - 25 -5.3 運行結(jié)果： - 26 -6 .M 次多項式曲線擬合 - 27 - # -6.1 算法說明 -27-6.2 算法程序： -28-7 . 不動點法解非線性方程 -33-7.1 算法說明

6、-33-7.3 運行結(jié)果： -35-8 . 二分法解非線性方程 -36-8.1 算法說明： -36-8.4 運行結(jié)果： -38-9，龍格 -庫塔法解微分方程 -39-9.1 算法說明： -39-9.2 算法程序： -40-9.3 運行界面： -41-10，遞歸梯形公式 -42-10.1 ，算法說明： -42-10.2 ，算法程序： -42-10.3 運行結(jié)果： -43-1.經(jīng)典四階龍格庫塔法解一階微分方程1.1 算法說明：龍格-庫塔(Runge-Kutta) 方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進行抑制，所以其實現(xiàn)原理也較復(fù)雜。該算法是構(gòu)建在數(shù)學(xué)支持的基

7、礎(chǔ)之上的。對于一精度的歐拉公式有：yi+1=yi+h*K1K1=f(xi,yi)當(dāng)用點 xi 處的斜率近似值K1 與右端點xi+1 處的斜率K2 的算術(shù)平均值作為平均斜率K* 的近似值，那么就會得到二階精度的改進歐拉公式：yi+1=yi+h*( K1+ K2)/2K1=f(xi,yi)K2=f(xi+h,yi+h*K1)依次類推，如果在區(qū)間xi,xi+1 內(nèi)多預(yù)估幾個點上的斜率值K1 、K2、Km,并用他們的加權(quán)平均數(shù)作為平均斜率K*的近似值，顯然能構(gòu)造出具有很高精度的高階計算公式。經(jīng)數(shù)學(xué)推導(dǎo)、求解， 可也就是在工程中應(yīng)用廣泛的經(jīng)典龍格數(shù)值計算課程設(shè)計庫塔算法：yi+1=yi+h*( K1+

8、2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)通常所說的龍格-庫塔法是指四階而言的，我們可以仿二階、三階的情形推導(dǎo)出常用的標(biāo)準(zhǔn)四階龍格-庫塔法公式1.2 ,算法程序：#include<iostream>#include<iomanip>#define M 10using namespace std;int main()long double feval(long double ,long double ,long double );lon

9、g double feval1(long double ,long double ,long double );long double f4,g4;long double h,a,b,xa,ya;long double xM+1,yM+1,TM+1;cout<<"請輸入?yún)^(qū)間左右端點a,b:"cin>>a>>b;cout<<endl;/xa,ya 是初值cout<<"請輸入函數(shù)在左端點的初值xa,ya:"cin>>xa>>ya;cout<<endl;/h是步長

10、h=(b-a)/M;x0=xa;y0=ya;for(int i=0;i<=M;i+)Ti=a+h*i; / 給 T 賦值for(int k=0;k<=M;k+)f0=feval(Tk,xk,yk);g0=feval1(Tk,xk,yk);f1=feval(Tk+h/2,xk+h/2*f0,yk+h/2*g0);g1=feval1(Tk+h/2,xk+h/2*f0,yk+h/2*g0);f2=feval(Tk+h/2,xk+h/2*f1,yk+h/2*g1);g2=feval1(Tk+h/2,xk+h/2*f1,yk+h/2*g1);f3=feval(Tk+h,xk+h*f2,yk+

11、h*g2);g3=feval1(Tk+h,xk+h*f2,yk+h*g2);xk+1=xk+h/6*(f0+2*f1+2*f2+f3);yk+1=yk+h/6*(g0+2*g1+2*g2+f3);cout<<setw(12)<<"T"<<setw(12)<<"x"<<setw(12)<<"y"<<endl;for(k=0;k<=M;k+)cout<<setw(12)<<Tk<<setw(12)<<x

12、k<<setw(12)<<yk<<endl;return 0;long double feval(long double t,long double x,long double y)long double f;f=x+2*y;return f;long double feval1(long double t,long double x,long double y)long double f;f=3*x+2*y;return f;- 11 -1.3,運行結(jié)果2.1 算法說明：2.高斯列主元法解線性方程組首先將線性方程組做成增廣矩陣，對增廣矩陣進行行變換。對于元

13、素ai ,i ，在第 i 列中，第i 行及以下的元素選取絕對值最大的元素，將該最大元素所在的行與第i 行交換 ， 然后采用高斯消元法用新得到的ai ,i 所在的第i 行消去第 i 行以下的元素。依次進行直到an,n 。從而得到上三角矩陣。2.2 算法程序：本程序包含enter.h ， remove.h， judge.h及 main.cpp 這四個函數(shù)1， main.cpp 函數(shù)#include<iostream.h>#include<math.h>#include "enter.h"/#include "remove.h"/#in

14、clude "judge.h"/int main()int row,col;cout<<" 請輸入方程的個數(shù):"cin>>row;cout<<" 請輸入未知數(shù)的個數(shù)加一:"cin>>col;double a100100,*p1,*p100;for(int i=0;i<row;i+)pi=ai;p1=p;enter(p1,row,col);/remove(p1,row,col);/int number=judge(p1,row,col); / if(number=0)包含名為ente

15、r 的頭函數(shù)包含名為remove的頭函數(shù)包含名為judge 的頭函數(shù)調(diào)用enter 函數(shù)調(diào)用remove函數(shù)調(diào)用judge 函數(shù)數(shù)值計算課程設(shè)計cout<<" 結(jié)果 "<<endl;cout<<" 此方程組無解"<<endl; / 無解時予以提示if(number=-1)cout<<" 結(jié)果 "<<endl;cout<<" 此方程組有無窮多解"<<endl; / 有無窮多解時予以提示if(number=1)cout<

16、;<" 結(jié)果 "<<endl;cout<<" 此方程組有唯一解"<<endl;/有唯一解時，解出方程組的解double x100;xcol-2=arow-1col-1/arow-1row-1; /回代法求解上三角矩陣for(i=col-3;i>=0;i-) double t=aicol-1;for(int j=col-2;j>i;j-) t=t-aij*xj;xi=t/aii;for(i=0;i<col-1;i+)cout<<"x"<<i+1<&

17、lt;"="<<xi<<"t" / 依次輸出方程組的解cout<<endl;return 0;2， enter.h 函數(shù)/ 輸入系數(shù)矩陣void enter(double *q,int m,int n)cout<<" 請按行輸入未知數(shù)前面的系數(shù)和等式右邊的常數(shù)："<<endl;for(int i=0;i<m;i+)for(int j=0;j<n;j+) cin>>qij;3， remove.h 函數(shù)/尋找列主元，并移動該行，最后化為上三角矩陣void

18、remove(double *q,int m,int n)int min=m;if(n-1<m)min=n-1;for(int i=0;i<min-1;i+)int k=i;double max=qii;for(int j=i+1;j<m;j+)if(fabs(qji)>fabs(max) max=qji;k=j;/找到第i 列從 aii 開始的絕對值最大的元素if(k!=i)for(int j=0;j<n;j+)double mat=qij;qij=qkj;qkj=mat;/通過換行以保證主對角線上的元素是其所在位置及以下元素中絕對值最大的一個int t=0;f

19、or(j=i;j<m;j+)if(qji=0)t+;/在消元前判斷aii 及其所在列以下元素是否都為零，不都為零再進行消元if(t!=m-i)for(int j=i+1;j<m;j+)for(int k=n-1;k>=i;k-) qjk=qjk-qik*qji/qii;for(i=0;i<n-2;i+)for(int j=i+1;j<m;j+)qji=0;/保證經(jīng)消元后得到上三角陣4， judge.h 函數(shù)/用增廣矩陣的秩判斷該矩陣所對應(yīng)的方程組的解的情況int judge(double *q,int m,int n)int r1=m;for(int i=0;i&

20、lt;m;i+)- 15 -數(shù)值計算課程設(shè)計-19 -int g=0;for(int j=0;j<n;j+)if(qij=0)g+;if(g=n)r1=r1-1;)int r2=m;for(i=0;i<m;i+)(int h=0;for(int j=0;j<n-1;j+)if(qij=0)h+;if(h=n-1)r2=r2-1;)if(r1=n-1&&r2=n-1)return 1;else if(r1=r2&&r1<n-1)return -1;else return 0;)/求出增廣矩陣的秩求出系數(shù)矩陣的秩據(jù)不同的秩的情況返回不同的值2

21、.3運行結(jié)果:3.牛頓法解非線性方程組3.1算法說明設(shè)R已知。第1步：計算函數(shù)fi(Pk,qk)f2(Pk,qk)P重復(fù)上述過程。F(Pk)第2步：計算雅可比矩陣J(R)第3步：求線性方程組J(Pk) P的解P。第4步：計算下一點Ri Rfi(Pk,qk) fi(Pk,qk) xyf2(Pk,qk) f2(Pk,qk) xyF(Pk)3.2算法程序:#include<iostream>#include<cmath>#define N 3/非線性方程組中方程個數(shù)、未知量個數(shù)#define Epsilon 0.0001 / 差向量1范數(shù)的上限#define Max 3/最大

22、迭代次數(shù)using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/計算向量函數(shù)的因變量向量 yyN陜西科技大學(xué)void ffjacobian(float xxN,float yyNN);/計算雅克比矩陣yyNNvoid inv_jacobian(float yyNN,float invNN);/計算雅克比矩陣的逆矩陣invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量x0 計算近似解向量x1float x0N

23、=0,0,0,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/ 如果取消對x0 的初始化，撤銷下面兩行的注釋符，就可以由鍵盤向x0 讀入初始近似解向量/for( i=0;i<N;i+)/cin>>x0i;cout<<" 初始近似解向量："<<endl;for (i=0;i<N;i+)cout<<x0i<<" " cout<<endl;cout<<endl;doiter=iter+1;cout

24、<<" 第 "<<iter<<" 次迭代開始"<<endl;/ 計算向量函數(shù)的因變量向量y0ff(x0,y0);/ 計算雅克比矩陣jacobianffjacobian(x0,jacobian);/ 計算雅克比矩陣的逆矩陣invjacobianinv_jacobian(jacobian,invjacobian);/ 由近似解向量x0 計算近似解向量x1newdundiedai(x0, invjacobian,y0,x1);/ 計算差向量的1 范數(shù) errornormerrornorm=0;for (i=0;i

25、<N;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornorm<Epsilon) break;for (i=0;i<N;i+) x0i=x1i; while (iter<Max);return 0;void ff(float xxN,float yyN)float x,y,z;int i;x=xx0;y=xx1;z=xx2;yy0=x*x-x+y*y+z*z-5;yy1=x*x+y*y-y+z*z-4;yy2=x*x+y*y+z*z+z-6;cout<<" 向量函數(shù)的因變量向量是："<&l

26、t;endl;for( i=0;i<N;i+)cout<<yyi<<" "cout<<endl;cout<<endl;void ffjacobian(float xxN,float yyNN)float x,y,z;int i,j;x=xx0;y=xx1;z=xx2;/jacobian have n*n elementyy00=2*x-1;yy01=2*y;yy02=2*z;yy10=2*x;yy11=2*y-1;yy12=2*z;yy20=2*x;yy21=2*y;yy22=2*z+1;cout<<&quo

27、t; 雅克比矩陣是："<<endl;for( i=0;i<N;i+)for(j=0;j<N;j+) cout<<yyij<<" "cout<<endl;cout<<endl;void inv_jacobian(float yyNN,float invNN)float augNN2,L;- # -陜西科技大學(xué)int i,j,k;cout<<"開始計算雅克比矩陣的逆矩陣："<<endl;for (i=0;i<N;i+) for(j=0;j<N

28、;j+)augij=yyij;for(j=N;j<N2;j+)if(j=i+N) augij=1;else augij=0;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=0;i<N;i+)for (k=i+1;k<N;k+)L=-augki/augii;for(j=i;j<N2;j+) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j

29、<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=N-1;i>0;i-)for (k=i-1;k>=0;k-)L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=

30、N-1;i>=0;i-)for(j=N2-1;j>=0;j-)augij=augij/augii;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout«augij«"cout«endl;for(j=N;j<N2;j+)invij-N=augij；)cout«endl;cout«"雅克比矩陣的逆矩陣："«endl;for (i=0;i<N;i+) for(j=0;j<N;j+)cout«invij«"cout

31、71;endl;)cout«endl;)void newdundiedai(float xON, float invNN,float yON,float x1N) (int ij;float sum=0;for(i=0;i<N;i+)sum=0;for(j=0;j<N;j+)sum=sum+invij*yOj;x1i=x0i-sum;)cout«"近似解向量："«endl;for (i=0;i<N;i+)cout«x1i«"cout«endl;cout«endl;)3.3運行結(jié)

32、果：-#-數(shù)值計算課程設(shè)計M /新M立堂關(guān)vt和逢試法De匕吧1牛迭代法用Ki結(jié)果為：x=-1.36628;y=-0.366281;z=2.36628初蛤Sffl晌量tJ M 0黯皤感%向量是:5T養(yǎng)克比矩陣是； -10 6 0 -1 G 0 0 1卜力鋁+算雅克比矩陣的逆矩陣: -1 QQ10a0 -1s01a0 Q 1 0 0 1-i Bsi0a0 -1q01ae n i s 0 1-1 0fi100(J Yg(d10-25 -數(shù)值計算課程設(shè)計4.龍貝格求積分算法4.1 算法說明生成J K的逼近表R(J,K),并以R(J 1,J 1)為最終解來逼近積分bf (x)dx R(J, J)a逼近

33、R(J,K)存在于一個特別的下三角矩陣中，第0列元素R(J,0)用基于2J個a,b子區(qū)間的連續(xù)梯形方法計算，然后利用龍貝格公式計算R(J,K)。當(dāng)1 K J時，第J行的元素為-31 -R(J,K) R(J,K 1)R(J,K 1) R(J 1,K 1)聲明在指定函數(shù)的指定點處的函數(shù)函數(shù)的循環(huán)當(dāng)|R(J,J) R(J 1,J 1)| tol時，程序在第(J 1)行結(jié)束。4.2 算法程序:/龍貝格方法求定積分/函數(shù)的修改可以在feval函數(shù)的定義中修改#include<iostream>#include<cmath>#include<iomanip> using

34、 namespace std;#define N 4 /定義全局變量 Nint main()long double feval(long double );/long double a,b,h;long double R1010=0;long double x;int M=1;long double err=1;int J=0;long double s=0;coutvv"請輸入兩個端點a,b:"cin>>a>>b;h=b-a;R00=h*(feval(a)+feval(b)while(err>0.0001)&&(J<N)

35、|(J<4) /J=J+1; h=h s=0;for(int p=1;p<=M;p+) (x=a+h*(2*p-1);s=s+feval(x);RJ0=RJ-10/2+h*s;M=2*M;for(int K=1;K<=J;K+) (RJK=RJK-1+(RJK-1-RJ-1K-1)/(4AK-1); err=fabs(RJ-1J-1-RJK); cout<<"該矩陣為："<<endl;for(int i=0;i<=J;i+)for(int j=0;j<=J;j+)(cout<<setiosflags(ios:

36、fixed)<<setprecision(8)<<setw(12)<<Rij; cout<<endl;cout<<"該式子在<<a<<","<<b<<"內(nèi)的積分為"<<RJJ;cout<<endl;return 0; long double feval(long double x) /函數(shù)的定義long double f;f=(x*x+x+1)*cos(x); return(f);4.3運行結(jié)果:5.三次緊壓樣條插

37、值5.1算法說明:三次緊壓樣條，確定SU), S(%)(如果導(dǎo)數(shù)已知，這是“最佳選擇”),m0(d0 S (x0)h0mN -一(S(Xn)5.2算法程序:#include<iostream>#include<iomanip>using namespace std;int main()const int n=4;const int N=n-1;/3long double Xn=0,1,2,3,Yn=0,0.5,2.0,1.5,dx0=0.2,dxn=-1,Hn-1,H1n-1,Dn-1,Un-2,AN-2,BN-1,CN-1,MN+1,temp,SN4;cout<

38、<"求解經(jīng)過點(0, 0.0),(1 , 0.5),(2 , 2.0)和(3, 1.5),而且一階導(dǎo)數(shù)邊界條件S'(0)=0.2和S'(3)=-1的三次壓緊樣條曲線："；for(int i=0;i<n-1;i+)Hi=Xi+1-Xi;for(int i1=0;i1<n-1;i1+)H1i1=Yi1+1-Yi1;for(int i2=0;i2<n-1;i2+)Di2=H1i2/Hi2;for(int i3=0;i3<n-1;i3+)Ui3=6*(Di3+1-Di3);for(int i4=0;i4<N-2;i4+)Ai4=H

39、i4+1;for(int i5=0;i5<N-1;i5+)Bi5=2*(Hi5+Hi5+1);for(int i6=0;i6<N-1;i6+)Ci6=Hi6+1;B0=B0-H0/2;U0=U0-3*(D0-dx0);BN-2=BN-2-HN-1/2;UN-2=UN-2-3*(dxn-DN-1);for(int k=1;k<N-1;k+)temp=Ak-1/Bk-1;Bk=Bk-temp*Ck-1;Uk=Uk-temp*Uk-1;MN-1=UN-2/BN-2;for(int k1=N-2;k1>0;k1-)Mk1=(Uk1-1-Ck1-1*Mk1+1)/Bk1-1;M0

40、=3*(D0-dx0)/H0-M1/2;MN=3*(dxn-DN-1)/HN-1-MN-1/2;for(int k2=0;k2<N;k2+)Sk20=(Mk2+1-Mk2)/(6*Hk2);Sk21=Mk2/2;Sk22=Dk2-Hk2*(2*Mk2+Mk2+1)/6;Sk23=Yk2;cout<<"S=:"<<endl;for(int j=0;j<N;j+)for(int j1=0;j1<4;j1+)cout<<setw(12)<<S皿1;cout<<endl;return 0;5.3運行結(jié)果:

41、數(shù)值計算課程設(shè)計6.M 次多項式曲線擬合6.1 算法說明設(shè)（xk,yk） N1有N個點，橫坐標(biāo)是確定的。最小二乘拋物線的系數(shù)表示為y f (x) Ax2 Bx C求解A, B和C的線性方程組為- 33 -Nxk4 A k1Nx3 Axk A k1Nx2 Axk A k1Nxk3 B k1Nxk2 B k1Nxk B k1NN22xk C ykxkk1k1NNxk Cykxkk1k1NNC ykk1數(shù)值計算課程設(shè)計6.2算法程序:#include<iostream>#include<cmath>#include<iomanip>using namespace

42、std;const int n=6,m=5,N=6,N2=2*N;int main()long double inv_jacobian(long double yyNN,10ng double invNN); long double Xn=0.25,1.0,1.5,2.0,2.4,5.0,Yn=23.1,1.68,1.0,0.84,0.826,1.2576,Bm+1,Fnm+1=0,F1m+1n,F2nm+1,Am+1m+1,Y1n,invm+1m+1,Cm+1,B1m+1;for(int i=0;i<m+1;i+)for(int j=0;j<n;j+)Fji=pow(Xj,i);

43、/F的轉(zhuǎn)置for(int i1=0;i1<m+1;i1+)for(int j1=0;j1<n;j1+)F1i1j1=Fj1i1;/輸出F1cout<<"F 的轉(zhuǎn)置 F'=:"<<endl;for(int p1=0;p1<m+1;p1+)for(int k1=0;k1<n;k1+)cout<<setw(13)<<F1p1k1;cout<<endl;cout<<endl;/輸出Fcout<<"F=:"<<endl;for(int p

44、2=0;p2<n;p2+)for(int k2=0;k2<m+1;k2+)cout<<setw(13)<<Fp2k2;cout<<endl;cout<<endl;cout<<"<以下計算矩陣A>"<<endl;/F的轉(zhuǎn)置F1乘以F得A矩陣for(int i2=0;i2<m+1;i2+) for(int j2=0;j2<n;j2+)for(int k1=0;k1<m+1;k1+)F2j2k1=F1i2j2*Fj2k1;/每次輸出F2cout<<&quo

45、t;第"<<i2+1<<"次乘法："<<endl;for(int p3=0;p3<n;p3+)for(int k3=0;k3<m+1;k3+)cout<<setw(13)<<F2p3k3;cout<<endl; cout<<endl;/得到矩陣Afor(int j3=0;j3<m+1;j3+) Ai2j3=0;/必須給此數(shù)賦初值for(int i3=0;i3<n;i3+) Ai2j3=Ai2j3+F2i3j3;/輸出Acout<<"A=

46、F'*F=:"<<endl;for(int p=0;p<m+1;p+)for(int k=0;k<m+1;k+)cout<<setw(13)<<Apk;cout<<endl;cout<<"<以下計算矩陣B>"<<endl;for(int i4=0;i4<m+1;i4+) for(int j4=0;j4<n;j4+)Y1j4=F1i4j4*Yj4;/每次輸出Y1cout<<"第"<<i4+1<<&

47、quot;次乘法："<<endl;for(int p4=0;p4<n;p4+)cout<<setw(13)<<Y1p4;cout<<endl;/得到矩陣BBi4=0;/必須給此數(shù)賦初值for(int i5=0;i5<m+1;i5+)Bi4=Bi4+Y1i5;cout<<"B=F'*Y'=:"<<endl;for(int i6=0;i6<n;i6+)cout<<setw(13)<<Bi6;cout<<endl;invm+1m+

48、1=inv_jacobian(A,inv);cout<<"A矩陣的逆矩陣："<<endl;for (i=0;i<N;i+) for(int j0=0;j0<N;j0+)cout<<setw(13)<<invij0;cout<<endl;cout<<"<以下計算矩陣C>"<<endl; for(int qq=0;qq<m+1;qq+) for(int w=0;w<n;w+)B1w=invqqw*Bw;/每次輸出B1cout<<

49、"第"<<qq+1<<"次乘法："<<endl;for(int ww=0;ww<n;ww+)cout<<setw(13)<<B1ww;cout<<endl;/得到矩陣CCqq=0;/必須給此數(shù)賦初值for(int e=0;e<m+1;e+)Cqq=Cqq+B1e;cout<<"C=A'*B="<<endl;for(int tt=0;tt<m+1;tt+)cout<<setw(13)<<Ctt

50、;cout<<endl;return 0;long double inv_jacobian(long double yyNN,10ng double invNN)long double augNN2,L;int i,j,k;cout<<"開始計算A矩陣的逆矩陣："<<endl;for (i=0;i<N;i+) for(j=0;j<N;j+)augij=yyij;for(j=N;j<N2;j+)if(j=i+N) augij=1;else augij=0;cout<<"第零次消去："<

51、<endl;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<"/"cout<<endl;cout<<endl;for (i=0;i<N;i+)for (k=i+1;k<N;k+)L=-augki/augii;for(j=i;j<N2;j+)augkj=augkj+L*augij;-39 -cout«"第一次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) co

52、utvvaugijvvr cout«endl;)cout«endl;for (i=N-1;i>0;i-)(for (k=i-1 ;k>=0;k-)L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij； )cout«"第二次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) coutvvaugijvvr cout«endl;)cout«endl;for (i=N-1;i>=0;i-) for

53、(j=N2-1;j>=0;j-)augiO=augij/augii;cout«"第三次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) coutvvaugijvvr cout«endl;for(j=N;j<N2;j+) invij-N=augij；) cout«endl;return(invm+1m+1);)6.3運行結(jié)果:n=F1 *F=：G12.153S.B725151.215bSQ.2443245.2212.159S.S72E151.SIG3245-221E892.E

54、38.0725151.215B8B.2443245.2215892.5787297151.215680.2443245.22158?2_573?2? . 739ZB98&80.2443245.2215892.57872?.73930081.9S632e+的63245.221E8T2.578729,7393Q81.95632e*00t?.?7385e*S06pl®懾矩酗23 .11.6g10.6261.2576酬欣乘法：5.7/1.661_51.60l.?8246,260情啾乘法：1-443751.6U2.25J.364-VS7763L.44用欷乘法： 0.36U93S1.6S3.3