【正文】 為 :endl。im。 //system(pause)。 } void load() {int i,j。im。jm+1。}數(shù)值計算課程設(shè)計 9 1122( , ) ( , )()( , ) ( , )k k k kkk k k kf p q f p qxyJPf p q f p qxy??????? ????12( , )() ( , )kkk kkf p qFP f p q??? ????1kkP P P? ? ?? 牛頓法解非線性方程組 、 算法說明 設(shè) kP 已知。 第 4 步：計算下一點 （ 33） 重復(fù)上述過程。 圖 32 牛頓法解非線性方程組算法程序運行結(jié)果 數(shù)值計算課程設(shè)計 11 、 牛頓法解非線性方程組算法程序代碼 includeiostream includecmath define N 2 // 非線性方程組中方程個數(shù)、未知量個數(shù) define epsilon // 差向量 1 范數(shù)的上限 define max 10 //最大迭代次數(shù) using namespace std。 int main() { void ff(float xx[N],float yy[N])。 void inv_jacobian(float yy[N][N],float inv[N][N])。 float x0[N]={,},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm。 cout初始解向量： endl。iN。 coutendl。 cout第 iter 次迭代開始： endl。 //jis uan jacobian ju zhen jacobian ffjacobian(x0,jacobian)。 //you jie xiang liang x0 ji suan jie xiang liang x1 newdim(x0, invjacobian,y0,x1)。 for (i=0。i++) errornorm=errornorm+fabs(x1[i]x0[i])。 for (i=0。i++) x0[i]=x1[i]。 return 0。 int i。 y=xx[1]。 yy[1]=x*x+4*y*y4。 for( i=0。i++) coutyy[i] 。 coutendl。 int i,j。 y=xx[1]。 yy[0][1]=1。 yy[1][1]=8*y。 for( i=0。i++) {for(j=0。j++) coutyy[i][j] 。 } coutendl。 int i,j,k。 for (i=0。i++) {for(j=0。j++) aug[i][j]=yy[i][j]。jN2。 else aug[i][j]=0。iN。jN2。 coutendl。 for (i=0。i++) { for (k=i+1。k++) {L=aug[k][i]/aug[i][i]。jN2。 } } for (i=0。i++) {for(j=0。j++) coutaug[i][j] 。 } coutendl。i0。k=0。 for(j=N21。j) aug[k][j]=aug[k][j]+L*aug[i][j]。iN。jN2。 coutendl。 for (i=N1。i) for(j=N21。j) aug[i][j]=aug[i][j]/aug[i][i]。iN。jN2。 牛頓法解非線性方程組 14 coutendl。jN2。 } coutendl。 for (i=0。i++) { for(j=0。j++) coutinv[i][j] 。 } coutendl。 //計算非線性方程組的近似解向量 float sum=0。iN。 for(j=0。j++) sum=sum+inv[i][j]*y0[j]。 } cout近似解向量： endl。iN。 } 數(shù)值計算課程設(shè)計 15 龍貝格求積分算法 、 算法說明 生成 JK? 的逼近表 ( , )RJ K ，并以 ( 1, 1)R J J??為最終解來逼近積分 ( ) ( , )ba f x dx R J J?? （ 41） 逼近 ( , )RJ K 存在于一個特別的下三角矩陣中，第 0列元素 ( ,0)RJ 用基于 2J 個 [a,b]子區(qū)間的連續(xù)梯形方法計算，然后利用龍貝格公式計算 ( , )RJ K 。 、 龍貝格求積分算法流程圖 圖 41 算法流程圖 龍貝格求積分算法 16 、 龍貝格求積分算法程序調(diào)試 我們以求解積分 ? ? )sin( 2xxf ? ，精度為 ，最高迭代 10 次為例，對所編寫的龍貝格求積分算法程序進行編譯和鏈接，經(jīng)執(zhí)行后得如下所示的窗口 圖 42龍貝格求 積分算法程序調(diào)試 說明：應(yīng)用 Romberg 算法求 ? ? )sin( 2xxf ? 在 ? ?21 區(qū)間上的精度為 的積分為。 define f(x) sin(x*x) //舉例函數(shù) define epsilon //精度 define MAXREPT 10 //迭代次數(shù) ,到最后仍達不到精度要求 ,則輸出 T(m=10). double Romberg(double aa, double bb) { //aa,bb 積分上下限 int m, n。 double s, q。 //精度要求 double *y = new double[MAXREPT]。//p 總是指示待計算元素的前一個元素 (同一行 ) //迭代初值 數(shù)值計算課程設(shè)計 17 h = bb aa。 m = 1。 ep = epsilon + 。amp。 for (int i=0。 i++)//求 Hn { x = aa + (i+)*h。 } p = (y[0] + h*p)/。 for (int k=1。 k++) { s = *s。 y[k1] = p。 龍貝格求積分算法 18 } p = fabs(q y[m1])。 y[m1] = q。 h = h/。 } int main() { double a,b。 cinab。 return 0。 y1=*(x1 1).^3 + *(x1 0).^3 *(x1 1) + *(x1 0)。 y2=*(x2 2).^3 *(x2 1).^3 *(x2 2) + *(x2 1)。 y3=*(x3 3).^3 + *(x3 2).^3 *(x3 3) + *(x3 2)。 Y=[0 2 ]。*39。S139。S239。S339。 const int max = 50。 float c[max], a[max], fxym[max]。 float b = (y[x2] y[x1]) / (x[x2] x[x1])。 } //求差分 void cal_m(int n) { //用追趕法求解出彎矩向量 M?? float B[max]。 for(int i = 1。 i++) B[i] = c[i] / (2 a[i]*B[i1])。 for(i = 1。 i++) fxym[i] = (fxym[i] a[i]*fxym[i1]) / (2 a[i]*B[i1])。 i = 0。 三次樣條插值算法 22 } void printout(int n)。 char ch。 cinn。 i = n。:39。 coutPlease put in Yi39。 ciny[i]。 i n。 coutPlease 輸入邊界條件 \n 1 :已知兩端的一階導(dǎo)數(shù) \n 2 :兩端的二階導(dǎo)數(shù)已知 \n 默認 :自然邊界條件 \n。 float f0, f1。 switch(t) { case 1:cout輸入 Y0\39。\n。 數(shù)值計算課程設(shè)計 23 c[0] = 1。 fxym[0] = 6*((y[1] y[0]) / (x[1] x[0]) f0) / h[0]。 break。 cinf0f1。 fxym[0] = 2*f0。 break。//待定 }。 i n。 for(i = 1。 i++) { a[i] = h[i1] / (h[i] + h[i1])。 }