北航數(shù)值分析大作業(yè)第三題(fortran)

1、“數(shù)值分析“計(jì)算實(shí)習(xí)大作業(yè)第三題SY1415215孔維鵬一、計(jì)算說明1、將 ， 分別代入方程組（A.3）得到關(guān)于t，u，v，w的的方程組，調(diào)用離散牛頓迭代子函數(shù)求出與 ， 對應(yīng)的 ， 。2、調(diào)用分片二次代數(shù)插值子函數(shù)在點(diǎn)( 處插值得到 ，得到數(shù)表 。3、對于 ，分別調(diào)用最小二乘擬合子函數(shù)計(jì)算系數(shù)矩陣 及誤差 ，直到滿足精度，即求得最小的k值及系數(shù)矩陣 。4、將 ， 分別代入方程組（A.3）得到關(guān)于 ， ， ， 的的方程組，調(diào)用離散牛頓迭代子函數(shù)求出與 ， 對應(yīng)的 ， ，調(diào)用分片二次代數(shù)插值子函數(shù)在點(diǎn)( 處插值得到 ；調(diào)用步驟3中求得的系數(shù)矩陣 求得 ，打印數(shù)表 。二、源程序（FORTRAN）P

2、ROGRAM SY1415215DIMENSION X(11,Y(21,T(6,U(6,Z(6,6,UX(11,21,TY(11,21,FXY(11,21,C(6,6DIMENSION X1(8,Y1(5,FXY1(8,5,PXY1(8,5,UX1(8,5,TY1(8,5REAL(8 X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE=第三題計(jì)算結(jié)果.TXTDO I=1,11X(I=0.08*(I-1ENDDODO I=1,21Y(I=0.5+0.05*(I-1ENDDO!*求解非線性方程組，得到z=f（t，u）的函數(shù)*DO

3、I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I,Y(J,UX(I,J,TY(I,JENDDOENDDO !*分片二次插值得到z=f（x，y）*DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J,TY(I,J,FXY(I,JENDDOENDDOWRITE (1,(數(shù)表(x,y,f(x,y:WRITE (1,(3X,X,7X,Y,10X,F(X,YDO I=1,11DO J=1,21WRITE(1,(1X,F5.2,2X,F5.3,2X,E20.13 X(I,Y(J,FXY(I,JENDDOWRITE (1,(ENDDO!

4、*最小二乘擬合得到P（x，y）*N=11M=21WRITE (1,( ,K和分別為：DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,EWRITE (1,(I3,2X,E20.13 K-1,EIF(E10E-7 EXITENDDOWRITE(1,( WRITE(1,(系數(shù)矩陣Crs（按行）為：DO I=1,KDO J=1,KWRITE (1,(E20.13,2X, C(I,JENDDOWRITE (1,(WRITE (*,(ENDDODO I=1,8X1(I=0.1*IENDDODO J=1,5Y1(J=0.5+0.2*JENDDODO I=1,8DO J=1

5、,5CALL DISNEWTON_NONLINEAR(X1(I,Y1(J,UX1(I,J,TY1(I,JENDDOENDDODO I=1,8DO J=1,5CALL INTERPOLATION(UX1(I,J,TY1(I,J,FXY1(I,JENDDOENDDOPXY1=0DO I=1,8DO J=1,5DO II=1,KDO JJ=1,KPXY1(I,J=PXY1(I,J+C(II,JJ*(X1(I*(II-1*(Y1(J*(JJ-1ENDDOENDDOENDDOENDDOWRITE(1,( WRITE(1,(數(shù)表（x,y,f(x,y,p(x,y:WRITE(1,(2X,X,6X,Y,12X

6、,F(X,Y,14X,P(X,YDO I=1,8DO J=1,5WRITE(1,(F5.3,2X,F5.3,2X,E20.13,2X,E20.13 X1(I,Y1(J,FXY1(I,J,PXY1(I,JENDDOWRITE (1,(ENDDOCLOSE (1END!*用離散牛頓法求解非線性方程組*SUBROUTINE DISNEWTON_NONLINEAR(X1,Y1,U,TPARAMETER (N=4REAL EPS !EPS為迭代精度，M為最大迭代次數(shù)DIMENSION X(N,H(N,Y(N,JA(N,N,E(N,XK(NREAL(8 JA,X,H,Y,E,XK,U,T,V,W,X1,Y

7、1,E1,E2F1(T,U,V,W=0.5*COS(T+U+V+W-X1-2.67F2(T,U,V,W=T+0.5*SIN(U+V+W-Y1-1.07F3(T,U,V,W=0.5*T+U+COS(V+W-X1-3.74F4(T,U,V,W=T+0.5*U+V+SIN(W-Y1-0.79EPS=10E-12M=100X=1.0DO K=1,MH=1 !計(jì)算Y=F（x）Y(1=F1(X(1,X(2,X(3,X(4Y(2=F2(X(1,X(2,X(3,X(4 Y(3=F3(X(1,X(2,X(3,X(4Y(4=F4(X(1,X(2,X(3,X(4!計(jì)算JA（N,N）E=0.0DO I=1,NDO J

8、=1,NDO JJ=1,NIF(JJ=J THENE(JJ=X(JJ+H(JJELSEE(JJ=X(JJENDIFENDDOIF(I=1 THENJA(I,J=(F1(E(1,E(2,E(3,E(4-F1(X(1,X(2,X(3,X(4/H(JELSEIF(I=2 THENJA(I,J=(F2(E(1,E(2,E(3,E(4-F2(X(1,X(2,X(3,X(4/H(JELSEIF(I=3 THENJA(I,J=(F3(E(1,E(2,E(3,E(4-F3(X(1,X(2,X(3,X(4/H(JELSEIF(I=4 THENJA(I,J=(F4(E(1,E(2,E(3,E(4-F4(X(1,X

9、(2,X(3,X(4/H(JENDIFENDDOENDDO !求解線性方程組CALL GAUSS(JA,XK,-Y,N!判斷精度CALL NORM(XK,N,E1 CALL NORM(X,N,E2IF(E1/E2TA.OR.(A(L,K=TA THENTA=A(L,KTL=LDO J=K,NT(K,J=A(K,JA(K,J=A(TL,JA(TL,J=T(K,JENDDOTB(K=B(KB(K=B(TLB(TL=TB(KENDIFENDDODO I=K+1,NM(I,K=A(I,K/A(K,KA(I,K=0DO J=K+1,NA(I,J=A(I,J-M(I,K*A(K,JENDDOB(I=B(I

10、-M(I,K*B(KENDDOENDDO!回代過程X(N=B(N/A(N,NDO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J*X(JENDDOX(K=(B(K-S/A(K,KENDDO RETURNEND!*求向量的無窮范數(shù)*SUBROUTINE NORM(X,N,ADIMENSION X(NREAL(8 X,AA=ABS(X(1DO I=2,NIF(ABS(X(IABS(X(I-1 THENA=ABS(X(IENDIFENDDORETURNEND!*分片二次代數(shù)插值*SUBROUTINE INTERPOLATION(U,V,W PARAMETER (N=6,M=6D

11、IMENSION X(N,Y(M,Z(M,N,LK(3,LR(3REAL(8 X,Y,Z,H,TREAL(8 U,V,W,LK,LR !U,V分別為插值點(diǎn)處的坐標(biāo)，W為插值結(jié)果INTEGER R!*數(shù)據(jù)賦值*DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&0.94,0.3,-0.18,-0.5,-0.66,-

12、0.66,&2.06,1.18,0.46,-0.1,-0.5,-0.74,&3.5,2.38,1.42,0.62,-0.02,-0.5/ H=0.4T=0.2!*計(jì)算K,R*IF(UX(N-1-H/2 THENK=N-1ELSEDO I=3,N-2IF(UX(I-H/2.AND.(U=X(I+H/2 THENK=IENDIFENDDOENDIFIF(VY(M-1-T/2 THENR=M-1ELSEDO J=3,M-2IF(VY(J-T/2.AND.(VN P=NIF(P20 P=20IF(QM Q=MIF(Q20 Q=20XX=0YY=0D1=NAPX(1=0.0DO I=1,NAPX(1=A

13、PX(1+X(IENDDOAPX(1=APX(1/D1DO J=1,MV(1,J=0.0DO I=1,NV(1,J=V(1,J+Z(I,JENDDOV(1,J=V(1,J/D1ENDDOIF(P1 THEND2=0.0APX(2=0.0DO I=1,NG=X(I-APX(1D2=D2+G*GAPX(2=APX(2+(X(I-XX*G*GENDDOAPX(2=APX(2/D2BX(2=D2/D1DO J=1,MV(2,J=0.0DO I=1,NG=X(I-APX(1V(2,J=V(2,J+Z(I,J*GENDDOV(2,J=V(2,J/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.

14、0APX(K=0.0DO J=1,MV(K,J=0.0ENDDODO I=1,NG1=1.0G2=X(I-APX(1DO J=3,KG=(X(I-APX(J-1*G2-BX(J-1*G1G1=G2G2=GENDDOD2=D2+G*GAPX(K=APX(K+X(I*G*GDO J=1,MV(K,J=V(K,J+Z(I,J*GENDDOENDDODO J=1,MV(K,J=V(K,J/D2ENDDOAPX(K=APX(K/D2BX(K=D2/D1D1=D2ENDDOD1=MAPY(1=0.0DO I=1,MAPY(1=APY(1+Y(IENDDOAPY(1=APY(1/D1DO J=1,PU(J,

15、1=0.0DO I=1,MU(J,1=U(J,1+V(J,IENDDOU(J,1=U(J,1/D1ENDDOIF(Q1THEND2=0.0APY(2=0.0DO I=1,MG=Y(I-APY(1D2=D2+G*GAPY(2=APY(2+(Y(I*G*GENDDOAPY(2=APY(2/D2BY(2=D2/D1DO J=1,PU(J,2=0.0DO I=1,MG=Y(I-APY(1U(J,2=U(J,2+V(J,I*GENDDOU(J,2=U(J,2/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K=0.0DO J=1,PU(J,K=0.0ENDDODO I=1,MG1=

16、1.0G2=Y(I-APY(1DO J=3,KG=(Y(I-APY(J-1*G2-BY(J-1*G1G1=G2G2=GENDDOD2=D2+G*GAPY(K=APY(K+Y(I*G*GDO J=1,PU(J,K=U(J,K+V(J,I*GENDDOENDDODO J=1,PU(J,K=U(J,K/D2ENDDOAPY(K=APY(K/D2BY(K=D2/D1D1=D2ENDDOV(1,1=1.0V(2,1=-APY(1V(2,2=1.0DO I=1,PDO J=1,QA(I,J=0.0ENDDOENDDODO I=3,QV(I,I=V(I-1,I-1V(I,I-1=-APY(I-1*V(I-1

17、,I-1+V(I-1,I-2IF(I=4 THENDO K=I-2,2,-1V(I,K=-APY(I-1*V(I-1,K+V(I-1,K-1-BY(I-1*V(I-2,KENDDOENDIFV(I,1=-APY(I-1*V(I-1,1-BY(I-1*V(I-2,1ENDDODO I=1,PIF(I=1 THENT(1=1.0T1(1=1.0ELSEIF(I=2 THENT(1=-APX(1T(2=1.0T2(1=T(1T2(2=T(2ELSET(I=T2(I-1T(I-1=-APX(I-1*T2(I-1+T2(I-2IF(I=4 THENDO K=I-2,2,-1T(K=-APX(I-1*T2

18、(K+T2(K-1-BX(I-1*T1(KENDDOENDIFT(1=-APX(I-1*T2(1-BX(I-1*T1(1T2(I=T(IDO K=I-1,1,-1T1(K=T2(KT2(K=T(KENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L=A(K,L+U(I,J*T(K*V(J,LENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(IDO J=1,MY1=Y(JX2=1.0DD=0.0DO K=1,PG=A(K,QDO KK=Q-1,1,-1G=G*Y1+A(K,KKENDDOG=G*X2DD=DD+GX2=X2*X1

19、ENDDODT=DD-Z(I,JDT1=DT1+DT*DTENDDOENDDORETURNEND三、計(jì)算結(jié)果數(shù)表(x,y,f(x,y:X Y UX TY F(X,Y0.00 0.500 1.345 0.243 0.4465040260217E+000.00 0.550 1.322 0.269 0.3246832410866E+000.00 0.700 1.255 0.350 0.3401912404030E-020.00 0.750 1.235 0.377 -0.8873582055387E-010.00 0.850 1.196 0.434 -0.2505346523772E+000.00 0

20、.900 1.177 0.463 -0.3202764787154E+000.00 0.950 1.159 0.492 -0.3826681228086E+000.00 1.000 1.142 0.521 -0.4377957495101E+000.00 1.050 1.125 0.550 -0.4857588795274E+000.00 1.100 1.109 0.580 -0.5266672522606E+000.00 1.150 1.093 0.609 -0.5606385181100E+000.00 1.200 1.079 0.639 -0.5877965434618E+000.00

21、1.250 1.064 0.669 -0.6082697644452E+000.00 1.300 1.050 0.699 -0.6221894327719E+000.00 1.350 1.037 0.729 -0.6296883895048E+000.00 1.400 1.024 0.759 -0.6308997691168E+000.00 1.450 1.011 0.790 -0.6259561598452E+000.08 0.500 1.415 0.228 0.6380152050067E+000.08 0.550 1.391 0.253 0.5066117295708E+000.08 0

22、.600 1.368 0.279 0.3821763157102E+000.08 0.650 1.346 0.306 0.2648635435247E+000.08 0.800 1.284 0.388 -0.4346806749373E-010.08 0.950 1.229 0.473 -0.2860046192736E+000.08 1.000 1.211 0.502 -0.3523860630008E+000.08 1.100 1.178 0.560 -0.4639049136278E+000.08 1.150 1.163 0.589 -0.5092367829793E+000.08 1.

23、200 1.148 0.619 -0.5477611307544E+000.08 1.250 1.133 0.649 -0.5795943699592E+000.08 1.300 1.119 0.679 -0.6048572233071E+000.08 1.350 1.106 0.709 -0.6236734340898E+000.08 1.400 1.093 0.739 -0.6361682760837E+000.08 1.450 1.080 0.769 -0.6424676727883E+000.08 1.500 1.068 0.799 -0.6426971176492E+000.16 0

24、.500 1.483 0.214 0.8400813766831E+000.16 0.550 1.460 0.239 0.6997641125864E+000.16 0.600 1.437 0.264 0.5660613612391E+000.16 0.650 1.414 0.290 0.4391716517306E+000.16 0.750 1.372 0.343 0.2063761909359E+000.16 0.850 1.333 0.398 0.2060692397577E-020.16 0.900 1.315 0.426 -0.8935398464083E-010.16 1.000

25、1.279 0.483 -0.2507999539520E+000.16 1.050 1.262 0.512 -0.3209321867511E+000.16 1.100 1.246 0.541 -0.3840977416274E+000.16 1.150 1.231 0.570 -0.4403822530309E+000.16 1.200 1.216 0.600 -0.4898811730659E+000.16 1.250 1.201 0.629 -0.5326979436866E+000.16 1.300 1.187 0.659 -0.5689418288271E+000.16 1.350

26、 1.174 0.689 -0.5987265475032E+000.16 1.400 1.161 0.718 -0.6221686739756E+000.16 1.450 1.148 0.748 -0.6393865589031E+000.16 1.500 1.136 0.778 -0.6504993634875E+000.24 0.550 1.527 0.225 0.9029273821703E+000.24 0.600 1.504 0.249 0.7605801844086E+000.24 0.650 1.482 0.275 0.6247152461364E+000.24 0.700 1

27、.460 0.301 0.4955197895447E+000.24 0.750 1.439 0.327 0.3731340505334E+000.24 0.800 1.419 0.354 0.2576567336524E+000.24 0.900 1.381 0.409 0.4764703563104E-010.24 0.950 1.363 0.437 -0.4684940702542E-010.24 1.100 1.313 0.523 -0.2885737112683E+000.24 1.150 1.297 0.552 -0.3554064608815E+000.24 1.250 1.26

28、7 0.610 -0.4689182248284E+000.24 1.350 1.240 0.669 -0.5561910731630E+000.24 1.400 1.227 0.699 -0.5902469746186E+000.24 1.500 1.202 0.759 -0.6397468483377E+000.32 0.600 1.570 0.236 0.9646076509368E+000.32 0.650 1.547 0.261 0.8203474277775E+000.32 0.700 1.526 0.286 0.6824477015560E+000.32 0.800 1.485

29、0.339 0.4263923704805E+000.32 0.850 1.466 0.365 0.3084629582599E+000.32 1.000 1.411 0.448 -0.4285967343901E-020.32 1.050 1.395 0.477 -0.9483383158541E-010.32 1.150 1.363 0.534 -0.2555538925125E+000.32 1.200 1.347 0.563 -0.3258401857129E+000.32 1.250 1.333 0.592 -0.3895069599390E+000.32 1.300 1.319 0

30、.621 -0.4466381259700E+000.32 1.350 1.305 0.650 -0.4973249497540E+000.32 1.450 1.279 0.710 -0.5797565124179E+000.40 0.700 1.591 0.272 0.8789600202794E+000.40 0.750 1.570 0.298 0.7391450727049E+000.40 0.800 1.550 0.324 0.6057448220594E+000.40 0.850 1.530 0.350 0.4788837887740E+000.40 0.900 1.512 0.37

31、7 0.3586506733933E+000.40 0.950 1.493 0.405 0.2451021277699E+000.40 1.100 1.443 0.488 -0.5525285639684E-010.40 1.200 1.412 0.545 -0.2220945000227E+000.40 1.250 1.397 0.574 -0.2956352193106E+000.40 1.300 1.383 0.603 -0.3626794362747E+000.40 1.350 1.369 0.632 -0.4233061773066E+000.40 1.400 1.356 0.662

32、 -0.4776011207937E+000.40 1.450 1.343 0.691 -0.5256554594343E+000.40 1.500 1.331 0.721 -0.5675647513312E+000.48 0.750 1.634 0.284 0.9363227839056E+000.48 0.800 1.613 0.310 0.7947044362728E+000.48 0.850 1.594 0.336 0.6593871570749E+000.48 0.900 1.575 0.363 0.5304876715431E+000.48 0.950 1.557 0.390 0.

33、4080885997266E+000.48 1.000 1.539 0.417 0.2922442321330E+000.48 1.100 1.506 0.472 0.8030848215007E-010.48 1.300 1.446 0.586 -0.2650070438339E+000.48 1.350 1.432 0.615 -0.3352378374822E+000.48 1.400 1.419 0.644 -0.3991645893867E+000.48 1.450 1.406 0.674 -0.4568681687555E+000.48 1.500 1.394 0.703 -0.5

34、084349841214E+000.56 0.800 1.676 0.297 0.9923495597011E+000.56 0.850 1.657 0.323 0.8490326592663E+000.56 1.000 1.602 0.402 0.4566585768956E+000.56 1.050 1.585 0.429 0.3386546704747E+000.56 1.100 1.568 0.457 0.2271082699420E+000.56 1.200 1.537 0.513 0.2339219534309E-010.56 1.250 1.522 0.541 -0.688186

35、3698047E-010.56 1.350 1.494 0.599 -0.2341526439382E+000.56 1.400 1.481 0.627 -0.3073911710726E+000.56 1.450 1.468 0.657 -0.3744348739671E+000.56 1.500 1.455 0.686 -0.4353605132598E+000.64 0.550 1.846 0.168 0.2034201225610E+010.64 0.900 1.699 0.336 0.9020461834743E+000.64 0.950 1.681 0.362 0.76322637

36、99936E+000.64 1.000 1.664 0.388 0.6306048550718E+000.64 1.050 1.646 0.415 0.5042529674728E+000.64 1.100 1.630 0.443 0.3842166978607E+000.64 1.150 1.614 0.470 0.2705203064766E+000.64 1.300 1.569 0.554 -0.3256648583588E-010.64 1.400 1.542 0.611 -0.2032515079566E+000.64 1.450 1.529 0.640 -0.27933038976

37、33E+000.64 1.500 1.516 0.669 -0.3493199222056E+000.72 0.500 1.931 0.139 0.2464684390094E+010.72 0.550 1.907 0.159 0.2278059075084E+010.72 0.600 1.884 0.181 0.2095251278336E+010.72 0.950 1.742 0.349 0.9536697141176E+000.72 1.050 1.707 0.402 0.6789308703755E+000.72 1.100 1.691 0.429 0.5507747946397E+0

38、00.72 1.150 1.675 0.456 0.4288254602727E+000.72 1.250 1.644 0.511 0.2036155380751E+000.72 1.350 1.616 0.568 0.3268547317272E-020.72 1.400 1.602 0.596 -0.8765320673369E-010.80 0.550 1.968 0.151 0.2526399604944E+010.80 0.600 1.945 0.172 0.2338411418941E+010.80 1.100 1.751 0.416 0.7259922777709E+000.80

39、 1.150 1.734 0.443 0.5961374446623E+000.80 1.200 1.719 0.470 0.4723864937793E+000.80 1.250 1.704 0.498 0.3547580832715E+000.80 1.300 1.689 0.525 0.2432542897086E+000.80 1.400 1.662 0.582 0.3855659402974E-010.80 1.450 1.649 0.610 -0.5469867387737E-01K和分別為：1 0.3220909652261E+012 0.4659944283092E-025 0.2547860927577E-07系數(shù)矩陣Crs（按行）為：-0.2773550687525E+00 -0.6923251420521E+00 0.9998496305197E+00 -0.7293194197118E+00 0.2642229227349E+00 -0.3850864968263E-01數(shù)表（x,y,f(x,y,p(x,y:X Y F(X,Y P(X,Y0.100 1.100 -0.4454976466566E+00 -0.4454989906626E+000.100 1.300 -0.5975666668168E