数值分析第三次大作业.docx

1、数值分析第三次大作业数值分析第三次大作业姓名：陈怡然学号：380911122010/6/16数值分析第三次大作业一、 设计方案：运用newton迭代法，求解已给方程的根对应的。并根据题中所给的u,t,z的二维数表对z进行插值得到，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，在所给的条件下停止运算。将带入z=f(x,y)和p(x,y)求出所求的值，并比较。二、源程序：#define N 6 /矩阵的阶；#define M 4 /方程组未知数个数； #define EPSL 1.0e-12 /迭代的精度水平；#define MAXDIGIT 1.0e-219 #define SIGSUM

2、 1.0e-7#define L 500 /迭代最大次数；#define K 10 /最小二乘法拟合时的次数上限； #define X_NUM 11#define Y_NUM 21#define OUTPUTMODE 1 /输出格式：0-输出至屏幕，1-输出至文件 #include #include double mat_uN = 0, 0.4, 0.8, 1.2, 1.6, 2.0, mat_tN = 0, 0.2, 0.4, 0.6, 0.8, 1.0, mat_ztuNN = -0.5, -0.34, 0.14, 0.94, 2.06, 3.5, -0.42, -0.5, -0.26,

3、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, , mat_val_zX_NUMY_NUM = 0, init_tuvw4 = 1,2,1,2, mat_crsKK = 0;FILE *p;int main()int i, j, x, y, kmin, tmp = 0;double tmpval;int g

4、etval_z(double, double, int, int);int leasquare();void result_out(int); if(OUTPUTMODE) p = fopen(c:Result.txt, w+); fprintf(p, 计算结果：n); fclose(p); for(i = 0; i X_NUM; i+) for(j = 0; j Y_NUM; j+) tmp+= getval_z(0.08 * i,0.5 + 0.05 * j, i, j); if(!tmp) printf(迭代求解 z=f(x,y) 完毕。n); else printf(迭代超过最大次数，

5、计算结果可能不准确！n); if(kmin=leasquare() printf(近似表达式已求出！n); for(i = 0; i 8; i+) for(j = 0; j 5; j+) tmp+= getval_z(0.1 * (i+1),0.5 + 0.2 * (j+1), i, j); if(!tmp) printf(迭代求解 z=f(x*,y*) 完毕。n); for(x = 0; x 8; x+) for(y = 0; y 5; y+) tmpval = 0; for(i = 0; i kmin; i+) for(j = 0; j kmin; j+) tmpval= tmpval +

6、 mat_crsij * pow(0.1 * (x+1), i) * pow(0.5 + 0.2 * (y+1), j); if(OUTPUTMODE) p = fopen(c:Result.txt, at+); fprintf(p, x*%d=%f, y*%d=%f: f(x*%d, y*%d)=%.12le, p(x*%d, y*%d)=%.12len, x+1, 0.1 * (x+1), y+1, 0.5 + 0.2 * (y+1), x+1, y+1, mat_val_zxy, x+1, y+1, tmpval); fclose(p); else printf(x%d=%f, y%d=

7、%f: f(x*%d, y*%d)=%.12le, p(x*%d, y*%d)=%.12len, x+1, 0.1 * (x+1), y+1, 0.5 + 0.2 * (y+1), x+1, y+1, mat_val_zxy, x+1, y+1, tmpval); printf(结果见C:Result.txt！); getchar(); else printf(近似表达式未求出，增大K值后重试n); getchar(); int getval_z(double x, double y, int i, int j)/牛顿迭代法求方程的解。double vec_tuvwM, vec_deltaM,

8、 fdaoMM, fM, mo_k, mo_delta_k, shang, t, u;int n = 0, k, flag_max, flag_stat; void solve(double fdaoMM, double vec_deltaM, double fM);double interpolation(double, double); flag_stat = 1; for(k = 0; k M; k+) vec_tuvwk = init_tuvwk; do f0 = -1.0 *(0.5 * cos(vec_tuvw0) + vec_tuvw1 + vec_tuvw2 + vec_tuv

9、w3 - x - 2.67); f1 = -1.0 *(vec_tuvw0 + 0.5 * sin(vec_tuvw1) + vec_tuvw2 + vec_tuvw3 - y - 1.07); f2 = -1.0 *(0.5 * vec_tuvw0 + vec_tuvw1 + cos(vec_tuvw2) + vec_tuvw3 - x - 3.74); f3 = -1.0 *(vec_tuvw0 + 0.5 * vec_tuvw1 + vec_tuvw2 + sin(vec_tuvw3) - y - 0.79); fdao00 = -0.5 * sin(vec_tuvw0); fdao01

10、 = 1; fdao02 = 1; fdao03 = 1; fdao10 = 1; fdao11 = 0.5 * cos(vec_tuvw1); fdao12 = 1; fdao13 = 1; fdao20 = 0.5; fdao21 = 1; fdao22 = -1 * sin(vec_tuvw2); fdao23 = 1; fdao30 = 1; fdao31 = 0.5; fdao32 = 1; fdao33 = cos(vec_tuvw3); solve(fdao, vec_delta, f); flag_max = 0; for(k = 1; k fabs(vec_deltaflag

11、_max) flag_max = k; mo_delta_k = vec_deltaflag_max; flag_max = 0; for(k = 1; k fabs(vec_tuvwflag_max) flag_max = k; mo_k = vec_tuvwflag_max; shang = fabs(mo_delta_k / mo_k); if(shang EPSL) t = vec_tuvw0 + vec_delta0, u = vec_tuvw1 + vec_delta1; flag_stat = 0; break; else for(k = 0; k M; k+) vec_tuvw

12、k+= vec_deltak; while(n+ = L); if(!flag_stat) if(OUTPUTMODE) p = fopen(c:Result.txt, at+); fprintf(p, x%d=%f, y%d=%f: f(x%d, y%d)=%.12len, i, 0.08*i, j, 0.5+0.05*j, i, j, mat_val_zij = interpolation(u,t); fclose(p); else printf(x%d=%f, y%d=%f: f(x%d, y%d)=%.12len, i, 0.08*i, j, 0.5+0.05*j, i, j, mat

13、_val_zij = interpolation(u,t); return flag_stat;void solve(double fdaoMM, double vec_deltaM, double fM)/消元法求方程的解int i, j, k, ik;double tmp, mik; for(k = 0; k M-1; k+) ik = k; for(i = k; i M; i+) if(fabs(fdaoikk) fabs(fdaoik) ik = i; for(j = k; j M; j+) tmp = fdaokj; fdaokj = fdaoikj; fdaoikj = tmp;

14、tmp = fk; fk = fik; fik = tmp; for(i = k + 1; i M; i+) mik = fdaoik / fdaokk; for(j = k; j = 0; k-) tmp = 0; for(j = k+1; j M; j+) tmp+= fdaokj * vec_deltaj; vec_deltak = (fk - tmp) / fdaokk; double interpolation(double u, double t)/二次插值int r, i, j, k;double coe_u3, coe_t3, z; z = 0; r = int(fabs(t

15、/ 0.2) + 0.5); k = int(fabs(u / 0.4) + 0.5); if(r = 0) r = 1; if(r = 5) r = 4; if(k = 0) k = 1; if(k = 5) k = 4; coe_u0 = (u - mat_uk)*(u - mat_uk + 1) / (mat_uk - 1- mat_uk) / (mat_uk - 1- mat_uk + 1); coe_u1 = (u - mat_uk - 1)*(u - mat_uk + 1) / (mat_uk- mat_uk - 1) / (mat_uk- mat_uk + 1); coe_u2

16、= (u - mat_uk - 1)*(u - mat_uk) / (mat_uk + 1- mat_uk - 1) / (mat_uk + 1- mat_uk); coe_t0 = (t - mat_tr)*(t - mat_tr + 1) / (mat_tr - 1- mat_tr) / (mat_tr - 1- mat_tr + 1); coe_t1 = (t - mat_tr - 1)*(t - mat_tr + 1) / (mat_tr- mat_tr - 1) / (mat_tr- mat_tr + 1); coe_t2 = (t - mat_tr - 1)*(t - mat_tr

17、) / (mat_tr + 1- mat_tr - 1) / (mat_tr + 1- mat_tr); for(i = r - 1; i = r + 1; i+) for(j = k - 1; j = k + 1; j+) z+= coe_ti - r + 1 * coe_uj - k + 1 * mat_ztuij; return z;int leasquare()/曲面拟合int x, y, i, j, k, k_max, k_now;double vec_xX_NUM, vec_yY_NUM, mat_btbKK = 0, mat_gtgKK = 0, mat_btbbtKX_NUM,

18、 mat_btbbtuKY_NUM, tmp, sigma;void inv(double matrixKK, int); k_now = 1; for(i = 0; i X_NUM; i+) vec_xi = 0.08 * i; for(j = 0; j Y_NUM; j+) vec_yj = 0.5 + 0.05 * j; do for(i = 0; i k_now; i+) for(j = 0; j k_now; j+) tmp = 0; for(k = 0; k X_NUM; k+) tmp+= pow(vec_xk, i) * pow(vec_xk, j); mat_btbij =

19、tmp; inv(mat_btb, k_now); for(i = 0; i k_now; i+) for(j = 0; j k_now; j+) tmp = 0; for(k = 0; k Y_NUM; k+) tmp+= pow(vec_yk, i) * pow(vec_yk, j); mat_gtgij = tmp; inv(mat_gtg, k_now); for(i = 0; i k_now; i+) for(j = 0; j X_NUM; j+) tmp = 0; for(k = 0; k k_now; k+) tmp+= mat_btbik * pow(vec_xj,k); ma

20、t_btbbtij = tmp; for(i = 0; i k_now; i+) for(j = 0; j Y_NUM; j+) tmp = 0; for(k = 0; k X_NUM; k+) tmp+= mat_btbbtik * mat_val_zkj; mat_btbbtuij = tmp; for(i = 0; i k_now; i+) for(j = 0; j k_now; j+) tmp = 0; for(k = 0; k Y_NUM; k+) tmp+= mat_btbbtuik * pow(vec_yk,j); mat_btbij = tmp; for(i = 0; i k_

21、now; i+) for(j = 0; j k_now; j+) tmp = 0; for(k = 0; k k_now; k+) tmp+= mat_btbik * mat_gtgkj; mat_crsij = tmp; sigma = 0; for(x = 0; x X_NUM; x+) for(y = 0; y Y_NUM; y+) tmp = 0; for(i = 0; i k_now; i+) for(j = 0; j k_now; j+) tmp+= mat_crsij * pow(0.08 * x, i) * pow(0.5 + 0.05 * y, j); sigma+= (tm

22、p - mat_val_zxy) * (tmp - mat_val_zxy); if(OUTPUTMODE) p = fopen(c:Result.txt, at+); fprintf(p, %d %.12lenn, k_now, sigma); fclose(p); else printf(%d %.12lenn, k_now, sigma); if(sigmaSIGSUM) if(OUTPUTMODE) p = fopen(c:Result.txt, at+); fprintf(p,Crs%d%d=nn, k_now, k_now); for(i = 0; i k_now; i+) fpr

23、intf(p,); for(j = 0; j k_now-1; j+) fprintf(p,%.12le, mat_crsij); fprintf(p,%.12le, mat_crsik_now-1); fprintf(p,n); fprintf(p,n); fclose(p); else printf(Crs%d%d=nn, k_now, k_now); for(i = 0; i k_now; i+) printf(); for(j = 0; j k_now-1; j+) printf(%.12le, mat_crsij); printf(%.12le, mat_crsik_now-1);

24、printf(,n); printf(nn); return k_now; while(k_now+ K); return 0; void inv(double matrixKK, int rank)/dolittle分解int i, j, k, t, n;double mat_tmpKK = 0, matrix_bakKK, vec_tmpK, vec_xK = 0, vec_dxK = 0, vec_rK = 0, tmp, max_x, max_dx;void preprocess(double KK, double K, int);void LUDecomposition(double

25、 KK, int);void LUSolve(double KK, double K, double K, int); n = 0; for(i = 0; i rank; i+) for(j = 0; j rank; j+) matrix_bakij = mat_tmpij = matrixij; vec_tmpi = 0; LUDecomposition(mat_tmp, rank); for(i = 0; i rank; i+) if(i=0) vec_tmpi = 1; else vec_tmpi - 1 = 0; vec_tmpi = 1; LUSolve(mat_tmp, vec_x, vec_tmp, rank); while(n+=3) for(j = 0; j rank; j+) tmp = 0; for(t = 0; t rank; t+) tmp+= matrix_bakjt * vec_xt; vec_rj = vec_tmpj - tmp; LUSolve(mat_tmp, vec_dx, vec_r, rank); max_x = vec_x0, max_dx = vec_dx0; for(j = 0; j