北航数值分析计算实习题目 第三题

/*Description:《数值分析》编程作业三 插值与逼近 实现及测试程序 *@date：2011-12-12 */ #include <iostream> #include <cmath> using namespace std; /* 分配矩阵空间 */ void allocateMemory(int m, int n, double **&p) { p = new double*[m]; for(int i=0; i<m; i++) { p[i] = new double[n]; memset(p[i], 0, sizeof(double)*n); } } /* 释放矩阵空间 */ void freeMemory(int m, double **&p) { for(int i=0; i<m; i++) delete []p[i]; delete []p; p = NULL; } /* 选出矩阵a的第k个列向量中，第k个分量到第n个分量绝对值最大分量所在的下标号 */ int selectMaxIk(double a[5][5], int n, int k) { double maxa; int maxk; maxa = fabs(a[k][k]); maxk = k; for(int i=k+1; i<=n; i++) { if(fabs(a[i][k]) > maxa) { maxa = fabs(a[i][k]); maxk = k; } } return maxk; } /* 交换矩阵a的第k行和第Ik行 */ void swapRows(double a[5][5], double b[5], int n, int k, int Ik) { double tmp; for(int j=k; j<=n; j++) //前k-1个元素全为0 { tmp = a[k][j]; a[k][j] = a[Ik][j]; a[Ik][j] = tmp; } tmp = b[k]; b[k] = b[Ik]; b[Ik] = tmp; } /* 列主元素高斯消去法求解方程组Ax=b */ void GaussSolver(double a[5][5], double b[5], int n, double x[5]) { //step 1，消元过程 double m; for(int k=1; k<=n-1; k++) { int Ik = selectMaxIk(a, n, k); if(Ik != k) swapRows(a, b, n, k, Ik); for(int i=k+1; i<=n; i++) { m = a[i][k]/a[k][k]; for(int j=k; j<=n; j++) { a[i][j] = a[i][j] - m*a[k][j]; } b[i] = b[i] - m*b[k]; } } //step 2，回代过程 double tmp = 0; x[n] = b[n]/a[n][n]; for(int k=n-1; k>=1; k--) { tmp = 0; for(int j=k+1; j<=n; j++) { tmp += a[k][j]*x[j]; } x[k] = (b[k] - tmp) / a[k][k]; } } /*计算向量x的无穷范数*/ double getVectorNorm(double x[5]) { double norm = fabs(x[1]); for(int i=2; i<=4; i++) if(norm < fabs(x[i])) norm = fabs(x[i]); return norm; } /* Newton迭代法求解非线性方程组F(x)=0 */ bool NewtonMethod(double result[5], double x0, double y0) { double Fderivative[5][5], Fminus[5]; double x[5] = {0, 1, 1, 1, 1}; //x的迭代初始值设为(1, 1, 1, 1) double detax[5] = {0}; const double precision = 1.0e-12; const double M = 10000; //迭代求解 for(int k=0; k<M;) { //计算F在x处的jocobi矩阵 Fderivative[1][1] = -0.5*sin(x[1]); Fderivative[1][2] = Fderivative[1][3] = Fderivative[1][4] = 1; Fderivative[2][2] = 0.5*cos(x[2]); Fderivative[2][1] = Fderivative[2][3] = Fderivative[2][4] = 1; Fderivative[3][3] = -sin(x[3]); Fderivative[3][1] = 0.5; Fderivative[3][2] = Fderivative[3][4] = 1; Fderivative[4][4] = cos(x[4]); Fderivative[4][2] = 0.5; Fderivative[4][1] = Fderivative[4][3] = 1; //计算F在x处的负向量 Fminus[1] = -(0.5*cos(x[1]) + x[2] + x[3] + x[4] - x0 - 2.67); Fminus[2] = -(x[1] + 0.5*sin(x[2]) + x[3] + x[4] - y0 - 1.07); Fminus[3] = -(0.5*x[1] + x[2] + cos(x[3]) + x[4] - x0 - 3.74); Fminus[4] = -(x[1] + 0.5*x[2] + x[3] + sin(x[4]) - y0 - 0.79); //用列主元素高斯消去法求解方程组Fderivative*detax=Fminus GaussSolver(Fderivative, Fminus, 4, detax); //是否已达精度要求？ if( getVectorNorm(detax)/getVectorNorm(x) <= precision) { //已达精度要求 for(int i=1; i<=4; i++) result[i] = x[i]; return true; } else { //还没有达到精度要求，继续迭代 for(int i=1; i<=4; i++) x[i] += detax[i]; k++; if(k == M) { cout << "Newton迭代法已达最大迭代次数，没有求出方程组的解。" << endl; return false; } }//if else }//for } /* 从z-ut数表中选出3*3个点进行分片二次代数差值，求得z=f(u,t)函数在某一点(u,t)处的值 */ double interpolation(double z_ut[6][6], double u0, double t0) { int i = (int)(u0/0.4 + 0.5); int j = (int)(t0/0.2 + 0.5); if(i==0) i = 1; if(i==5) i = 4; if(j==0) j = 1; if(j==5) j = 4; double u[6], t[6]; u[0] = t[0] = 0; for(int k=1; k<=5; k++) { u[k] = u[k-1] + 0.4; t[k] = t[k-1] + 0.2; } //插值公式参见课本P101公式（5.17） double sum = 0; for(int k=i-1; k<=i+1; k++) { for(int r=j-1; r<=j+1; r++) { double tmp1 = 1, tmp2 = 1; for(int w=i-1; w<=i+1; w++) { if(w==k) continue; tmp1 *= (u0 - u[w])/(u[k] - u[w]); } for(int w=j-1; w<=j+1; w++) { if(w==r) continue; tmp2 *= (t0 - t[w])/(t[r] - t[w]); } sum += tmp1 * tmp2 * z_ut[r][k]; } } return sum; } /* 确定方程z=f(x,y)在11*21个点对上的函数值 */ void getFuncfxy(double z_xy[11][21], int lowi, int highi, int lowj, int highj, double detax, double detay) { //z-ut数表 double z_ut[6][6] = {{-0.5, -0.34, 0.14, 0.94, 2.06, 3.5}, {-0.42, -0.5, -0.26, 0.3, 1.18, 2.38}, {-0.18, -0.5, -0.5, -0.18, 0.46, 1.42}, {0.22, -0.34, -0.58, -0.5, -0.1, 0.62}, {0.78, -0.02, -0.5, -0.66, -0.5, -0.02}, {1.5, 0.46, -0.26, -0.66, -0.74, -0.5}}; //对每一个点对(xi,yi)求函数值z=f(xi,yi) for(int i=lowi; i<=highi; i++) { for(int j=lowj; j<=highj; j++) { double xi = detax * i; double yj = 0.5 + detay * j; double xyMap2ut[5] = {0}; //Newton迭代法求解非线性方程组 NewtonMethod(xyMap2ut, xi, yj); //分片二次插值求解(ui,ti)处的z点值，也就是与（ui,ti)对应的点（xi,yi)的函数值z=f(xi,yi) z_xy[i][j] = interpolation(z_ut, xyMap2ut[2], xyMap2ut[1]); } } } /* 利用LU分解法对矩阵求逆，求得的逆矩阵存放回原来的矩阵中 */ void getReversedMatrix(double **matrix,int rank) { double **Lower,**Upper,**ReverseLower,**ReverseUpper; //分配空间，初始化临时矩阵 allocateMemory(rank, rank, Lower); allocateMemory(rank, rank, Upper); allocateMemory(rank, rank, ReverseLower); allocateMemory(rank, rank, ReverseUpper); for(int i=0;i<rank;i++) { for(int j=0;j<rank;j++) { Lower[i][j]=0.0; Upper[i][j]=0.0; ReverseLower[i][j]=0.0; ReverseUpper[i][j]=0.0; } } //开始进行LU分解 for(int i=0;i<rank;i++) Lower[i][i]=1.0; for(int k=0;k<rank;k++) { for(int j=k;j<rank;j++) { double sum=0.0; for(int t=0;t<=k-1;t++) { sum+=Lower[k][t]*Upper[t][j]; } Upper[k][j]=matrix[k][j]-sum; } for(int i=k+1;i<rank;i++) { double sum=0.0; for(int t=0;t<=k-1;t++) { sum+=Lower[i][t]*Upper[t][k]; } Lower[i][k]=matrix[i][k]-sum; Lower[i][k]=Lower[i][k]/Upper[k][k]; } } for(int i=0;i<rank;i++) { ReverseLower[i][i]=1/Lower[i][i]; ReverseUpper[i][i]=1/Upper[i][i]; } //下三角阵L求逆 for(int j=0;j<rank;j++) { for(int i=j+1;i<rank;i++) { double sum=0.0; for(int k=j;k<i;k++) { sum+=Lower[i][k]*ReverseLower[k][j]; } ReverseLower[i][j]=-sum/Lower[i][i]; } } //上三角阵U求逆 for(int j=rank-1;j>=1;j--) { for(int i=j-1;i>=0;i--) { double sum=0.0; for(int k=i+1;k<=j;k++) { sum+=Upper[i][k]*ReverseUpper[k][j]; } ReverseUpper[i][j]=-sum/Upper[i][i]; } } for(int i=0;i<rank;i++) { for(int j=0;j<rank;j++) { matrix[i][j]=0.0; for(int k=0;k<rank;k++) { matrix[i][j]+=ReverseUpper[i][k]*ReverseLower[k][j]; } } } //释放空间 freeMemory(rank, Lower); freeMemory(rank, Upper); freeMemory(rank, ReverseLower); freeMemory(rank, ReverseUpper); } /* 计算两个矩阵的乘积 */ void getMatrixProduct(int m, int n, int r, double **&A, double **&B, double **&C) { for(int i=0; i<m; i++) { for (int j=0; j<r; j++) { double tmp = 0; for(int k=0; k<n; k++) tmp += A[i][k]*B[k][j]; C[i][j] = tmp; } } } /* 计算系数矩阵C */ void getMatrixC(int n, double **&B, double **&G, double z_xy[11][21], double **&RBTB, double **&RGTG, double **&C) { double **RB, **U, **tmp1, **tmp2, **tmp3; allocateMemory(n, 11, RB); allocateMemory(11, 21, U); allocateMemory(n, 11, tmp1); allocateMemory(n, 21, tmp2); allocateMemory(n, n, tmp3); //对矩阵B转置 for(int i=0; i<11; i++) for(int j=0; j<n; j++) RB[j][i] = B[i][j]; for(int