【正文】 [n1]) / (x[n] x[n1])) / h[n1]。 cint。:39。 do{ cout輸入 x的最大下標(biāo) :。 i = n。 return (a b)/(x[x3] x[x1])。) gtext(39。 X=[0 1 2 3]。 cout積分結(jié)果 :Romberg(a, b)endl。 m = m + 1。//求 T2n = 1/2(Tn+Hn),用 p指示 //求第 m行元素 ,根據(jù) Romberg 計算表本行的前一 個元素 (p 指示 ), //和上一行左上角元素 (y[k1]指示 )求得 . s = 。 //迭代計算 while ((ep = epsilon) amp。 double ep。i++) coutx1[i] endl。i++) { sum=0。jN。 for(j=N。j=0。i++) {for(j=0。i) { for (k=i1。iN。iN。 } 數(shù)值計算課程設(shè)計 13 for (i=0。iN。jN。 yy[0][0]=2*x2。iN。 } 牛頓法解非線性方程組 12 void ff(float xx[N],float yy[N]) { float x,y。iN。 do { iter=iter+1。 void newdim(float x0[N], float inv[N][N],float y0[N],float x1[N])。 第 1 步：計算函數(shù) （ 31） 第 2 步：計算雅可比矩陣 （ 32） 第 3 步：求線性方程組 ( ) ( )kkJ P P F P? ? ? 的解 P? 。 return 0。jm。 for(r=i。 a[c][n]=s。im。 int c,k,n,p,r。 Km max j=i:max=|*(head+(m+1)*i+i)|。 對第元素 iia ，在第 i 列中，第 i行及以下的元素選取絕對值最大的元素，將該元素最大的行與第 i 行交換，然后采用高斯消元法將新得到的 iia 消去第 i行以下的元素。i=(ba)/step。 經(jīng)典四階龍格庫塔法解一階微分方程 4 //double x0=0,y0=1。} int main() {double f(double x, double y)。 x1=xnh。 double h=(xnx0)/step。由于此算法精度高，采取措施對誤差進行抑制，所以其實現(xiàn)原理也較復(fù)雜。 k2=f(x0+h/2, y0+h*k1/2)。 k2=f(x1+h/2, y1+h*k1/2)。 // int step。 cout請輸入步長： 。} return(0)。 再對得到的上三角矩陣進行回代操作，即可以得到方程組的解。 void load()。 cinm。j++) c=(fabs(a[j][i])fabs(a[i][i]))?j:i。p++) couta[i][p]\t。 } } x[m1]=a[m1][m]/a[m1][m1]。} cout該方程組的解為 :endl。im。 圖 32 牛頓法解非線性方程組算法程序運行結(jié)果 數(shù)值計算課程設(shè)計 11 、 牛頓法解非線性方程組算法程序代碼 includeiostream includecmath define N 2 // 非線性方程組中方程個數(shù)、未知量個數(shù) define epsilon // 差向量 1 范數(shù)的上限 define max 10 //最大迭代次數(shù) using namespace std。 cout初始解向量： endl。 //jis uan jacobian ju zhen jacobian ffjacobian(x0,jacobian)。 for (i=0。 y=xx[1]。 coutendl。 yy[1][1]=8*y。 } coutendl。j++) aug[i][j]=yy[i][j]。jN2。k++) {L=aug[k][i]/aug[i][i]。j++) coutaug[i][j] 。 for(j=N21。 coutendl。iN。 } coutendl。 } coutendl。j++) sum=sum+inv[i][j]*y0[j]。 、 龍貝格求積分算法流程圖 圖 41 算法流程圖 龍貝格求積分算法 16 、 龍貝格求積分算法程序調(diào)試 我們以求解積分 ? ? )sin( 2xxf ? ，精度為 ，最高迭代 10 次為例，對所編寫的龍貝格求積分算法程序進行編譯和鏈接，經(jīng)執(zhí)行后得如下所示的窗口 圖 42龍貝格求 積分算法程序調(diào)試 說明：應(yīng)用 Romberg 算法求 ? ? )sin( 2xxf ? 在 ? ?21 區(qū)間上的精度為 的積分為。//p 總是指示待計算元素的前一個元素 (同一行 ) //迭代初值 數(shù)值計算課程設(shè)計 17 h = bb aa。 for (int i=0。 k++) { s = *s。 h = h/。 y1=*(x1 1).^3 + *(x1 0).^3 *(x1 1) + *(x1 0)。*39。 const int max = 50。 for(int i = 1。 i = 0。 i = n。 i n。\n。 cinf0f1。 i n。 cout\n輸出三次樣條插值函數(shù)： \n。)。 if(t 0) coutt*(x x[i+1])^3。 else cout t*(x x[i+1])。 Y=[15 5 1 5]。擬合曲線 39。j++) if(i==j) E[i][j]=1。 } for(k=0。 } } } } int main() {int n,M,i,j,k。i++) cinX[i]。in。j++) { BF[j][i]=F[i][j]。 //計算 F 的轉(zhuǎn)置 BF 與 Y 的乘 for(i=0。i++) for(j=0。 cout\nP(x)=。 、 牛頓 拉弗森迭代解非線性方程 算法程序調(diào)試 圖 71牛頓 拉弗森迭代解非線性方程算法程序調(diào)試 、 牛頓 拉弗森迭代解非線性方程算法代碼 includeiostream using namespace std。k=max1。 cout方程的根的誤差估計為 :errendl。 int k=0,max1=0,i=0。k++) { P[k1]=g(P[k2])。 cout近似值的誤差為 :errendl。 double func(double X,int k,double x[],int n)。 int i。in。in。 define delta 設(shè)計體會及今后的改進意見 38 define n 3 define max1 1000 main() { int i,j,k。 coutPlease input matrix B(n*1)endl。 for (k=0。 X[j]=(X[j]+A[j][j]*P[j])/A[j][j]。in。 ” 千里之行始于足下 ” ，通過這次課程設(shè)計，我深深體會到這句千古名言的真正含義 。 由于 我 的設(shè)計能力有限，在設(shè)計過程中難免出現(xiàn)錯誤，懇請老師們多多指教 ,我十分樂意接受你們的批評與指正， 我 將萬分感謝。 在此感謝我們的 劉海峰 老師 .，老師嚴(yán)謹(jǐn)細致、一絲不茍的作風(fēng)一直是我工作、學(xué)習(xí)中的榜樣；老師循循善誘的教導(dǎo)和不拘一格的思路給予我無盡的啟迪；這次 數(shù)值分析課程設(shè)計 的每個實驗細節(jié)和每個數(shù)據(jù)，都離不開老師您的細心指導(dǎo)。 for (i=0。i++) norm=norm+pow((X[i]P[i]),2)。j++) { X[j]=B[j]。 coutPlease input matrix P(n*1)endl。in。 } return Pn。 } double func(double X,int k,double x[],int n) { int i。 cinx[i]y[i]。 cinn。i++) { coutsetw(10)P[i]。 if((errtol)||(relerrtol)) break。 max1=100。 } 數(shù)值計算課程設(shè)計 33 不動點法解非線性方程 、 算法說明 先將 ()fx改寫成 ()x gx? 然后對 ()x gx? 進行迭代，即 1()kkx g x ?? 其中 1,2,...k? 然后判斷 1||kkxx ????是否成立，成立則返回 kx ，不成立就重復(fù)以上步驟 、 不動點法解非線性方程算法程序調(diào)試 我們將編寫好的不動點迭代法解非線性方程算法程序進行編譯，鏈接和執(zhí)行后得如下所示結(jié)果。 p0=p1。 cout牛頓拉弗森法解非線性方程 f(x)=x^22x1endl。 else cout+C[i]*x^i。 for(i