数值分析第三次大作业

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,

26、 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 max_x) max_x = vec_xj; if(vec_dxj max_dx) max_dx = vec_dxj; for(j = 0; j ran

27、k; j+) vec_xj+= vec_dxj; for(j = 0; j rank; j+) matrixji = vec_xj; void preprocess(double matrixkk, double vec_tmpk, int rank)/改变条件数int i, k;double tmp; for(i = 0; i rank; i+) tmp = matrixi0; for(k = 0; k tmp) tmp = matrixik; for(k = 0; k rank; k+) matrixik/= tmp; vec_tmpi/= tmp; void ludecompositio

28、n(double matrixkk, int rank)/lu分解int i, j, k, t;double tmp; for(k = 0; k rank; k+) for(i = k; i rank; i+) tmp = 0; for(t = 0; t k; t+) tmp+= matrixit * matrixtk; matrixik-= tmp; if(k rank - 1) for(j = k + 1; j rank; j+) tmp = 0; for(t = 0; t k; t+) tmp+= matrixkt * matrixtj; matrixkj = (matrixkj - t

29、mp) / matrixkk; else void lusolve(double matrixkk, double vec_xk, double vec_tmpk, int rank)/求解矩阵的逆int i, j, t;double tmp; vec_x0 = vec_tmp0 / matrix00; for(i = 1; i rank; i+) tmp = 0; for(t = 0; t = 0; i-) tmp = 0; for(t = i + 1; t rank; t+) tmp+= matrixit * vec_xt; vec_xi-= tmp; 计算结果：x0=0.000000,

30、y0=0.500000: f(x0, y0)=4.465040184799e-01x0=0.000000, y1=0.550000: f(x0, y1)=3.246832629274e-01x0=0.000000, y2=0.600000: f(x0, y2)=2.101596866825e-01x0=0.000000, y3=0.650000: f(x0, y3)=1.030436083159e-01x0=0.000000, y4=0.700000: f(x0, y4)=3.401895562659e-03x0=0.000000, y5=0.750000: f(x0, y5)=-8.8735

31、81363801e-02x0=0.000000, y6=0.800000: f(x0, y6)=-1.733716327497e-01x0=0.000000, y7=0.850000: f(x0, y7)=-2.505346114666e-01x0=0.000000, y8=0.900000: f(x0, y8)=-3.202765063876e-01x0=0.000000, y9=0.950000: f(x0, y9)=-3.826680661097e-01x0=0.000000, y10=1.000000: f(x0, y10)=-4.377957667384e-01x0=0.000000

32、, y11=1.050000: f(x0, y11)=-4.857589414438e-01x0=0.000000, y12=1.100000: f(x0, y12)=-5.266672548835e-01x0=0.000000, y13=1.150000: f(x0, y13)=-5.606384797965e-01x0=0.000000, y14=1.200000: f(x0, y14)=-5.877965387677e-01x0=0.000000, y15=1.250000: f(x0, y15)=-6.082697790899e-01x0=0.000000, y16=1.300000:

33、 f(x0, y16)=-6.221894528764e-01x0=0.000000, y17=1.350000: f(x0, y17)=-6.296883781856e-01x0=0.000000, y18=1.400000: f(x0, y18)=-6.308997600028e-01x0=0.000000, y19=1.450000: f(x0, y19)=-6.259561525454e-01x0=0.000000, y20=1.500000: f(x0, y20)=-6.149885466094e-01x1=0.080000, y0=0.500000: f(x1, y0)=6.380

34、152265102e-01x1=0.080000, y1=0.550000: f(x1, y1)=5.066117551462e-01x1=0.080000, y2=0.600000: f(x1, y2)=3.821763692772e-01x1=0.080000, y3=0.650000: f(x1, y3)=2.648634911536e-01x1=0.080000, y4=0.700000: f(x1, y4)=1.547802002848e-01x1=0.080000, y5=0.750000: f(x1, y5)=5.199268349093e-02x1=0.080000, y6=0

35、.800000: f(x1, y6)=-4.346804020491e-02x1=0.080000, y7=0.850000: f(x1, y7)=-1.316010567885e-01x1=0.080000, y8=0.900000: f(x1, y8)=-2.124310883088e-01x1=0.080000, y9=0.950000: f(x1, y9)=-2.860045510580e-01x1=0.080000, y10=1.000000: f(x1, y10)=-3.523860789794e-01x1=0.080000, y11=1.050000: f(x1, y11)=-4

36、.116554565222e-01x1=0.080000, y12=1.100000: f(x1, y12)=-4.639049115188e-01x1=0.080000, y13=1.150000: f(x1, y13)=-5.092367247005e-01x1=0.080000, y14=1.200000: f(x1, y14)=-5.477611179623e-01x1=0.080000, y15=1.250000: f(x1, y15)=-5.795943883391e-01x1=0.080000, y16=1.300000: f(x1, y16)=-6.048572588895e-

37、01x1=0.080000, y17=1.350000: f(x1, y17)=-6.236734213318e-01x1=0.080000, y18=1.400000: f(x1, y18)=-6.361682484133e-01x1=0.080000, y19=1.450000: f(x1, y19)=-6.424676566900e-01x1=0.080000, y20=1.500000: f(x1, y20)=-6.426971026996e-01x2=0.160000, y0=0.500000: f(x2, y0)=8.400813957651e-01x2=0.160000, y1=

38、0.550000: f(x2, y1)=6.997641656726e-01x2=0.160000, y2=0.600000: f(x2, y2)=5.660614423514e-01x2=0.160000, y3=0.650000: f(x2, y3)=4.391716081175e-01x2=0.160000, y4=0.700000: f(x2, y4)=3.192421380408e-01x2=0.160000, y5=0.750000: f(x2, y5)=2.063761923874e-01x2=0.160000, y6=0.800000: f(x2, y6)=1.006385238914e-01x2=0.160000, y7=0.850000: f(x2, y7)=2.060740067835e-03x2=0.160000, y8=0.900000: f(x2, y8)=-8.93