北航数值分析报告大作业第三题fortran.docx

1、北航数值分析报告大作业第三题fortran“数值分析“计算实习大作业第三题SY1415215孔维鹏一、计算说明1、将，分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与，对应的，。2、调用分片二次代数插值子函数在点()处插值得到)，得到数表。3、对于，分别调用最小二乘拟合子函数计算系数矩阵及误差，直到满足精度，即求得最小的k值及系数矩阵。4、将，分别代入方程组（A.3）得到关于，的的方程组，调用离散牛顿迭代子函数求出与，对应的，调用分片二次代数插值子函数在点()处插值得到)；调用步骤3中求得的系数矩阵求得，打印数表。二、源程序（FORTRAN）PROGRAM

2、 SY1415215DIMENSION X(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21),C(6,6)DIMENSION X1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE=第三题计算结果.TXT) DO I=1,11 X(I)=0.08*(I-1)ENDDODO I=1,21 Y(I)=0.5+0.05*(I-1)ENDDO!*求解非线性方程组

3、，得到z=f（t，u）的函数*DO I=1,11 DO J=1,21 CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J) ENDDOENDDO !*分片二次插值得到z=f（x，y）*DO I=1,11 DO J=1,21 CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J) ENDDOENDDOWRITE (1,(数表(x,y,f(x,y):)WRITE (1,(3X,X,7X,Y,10X,F(X,Y)DO I=1,11 DO J=1,21 WRITE(1,(1X,F5.2,2X,F5.3,2X,E20.13)

4、 X(I),Y(J),FXY(I,J) ENDDO WRITE (1,()ENDDO!*最小二乘拟合得到P（x，y）*N=11M=21WRITE (1,( ,K和分别为：)DO K=1,20 CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E) WRITE (1,(I3,2X,E20.13) K-1,E IF(E10E-7) EXITENDDOWRITE(1,( )WRITE(1,(系数矩阵Crs（按行）为：)DO I=1,K DO J=1,K WRITE (1,(E20.13,2X,) C(I,J) ENDDO WRITE (1,() WRITE (*,()ENDDODO

5、I=1,8 X1(I)=0.1*IENDDODO J=1,5 Y1(J)=0.5+0.2*JENDDODO I=1,8 DO J=1,5 CALL DISNEWTON_NONLINEAR(X1(I),Y1(J),UX1(I,J),TY1(I,J) ENDDOENDDODO I=1,8 DO J=1,5 CALL INTERPOLATION(UX1(I,J),TY1(I,J),FXY1(I,J) ENDDOENDDOPXY1=0DO I=1,8 DO J=1,5 DO II=1,K DO JJ=1,K PXY1(I,J)=PXY1(I,J)+C(II,JJ)*(X1(I)*(II-1)*(Y1(

6、J)*(JJ-1) ENDDO ENDDO ENDDOENDDOWRITE(1,( )WRITE(1,(数表（x,y,f(x,y),p(x,y):)WRITE(1,(2X,X,6X,Y,12X,F(X,Y),14X,P(X,Y)DO I=1,8 DO J=1,5 WRITE(1,(F5.3,2X,F5.3,2X,E20.13,2X,E20.13) X1(I),Y1(J),FXY1(I,J),PXY1(I,J) ENDDO WRITE (1,()ENDDO CLOSE (1)END!*用离散牛顿法求解非线性方程组*SUBROUTINE DISNEWTON_NONLINEAR(X1,Y1,U,T)

7、PARAMETER (N=4)REAL EPS !EPS为迭代精度，M为最大迭代次数DIMENSION X(N),H(N),Y(N),JA(N,N),E(N),XK(N)REAL(8) JA,X,H,Y,E,XK,U,T,V,W,X1,Y1,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,M H

8、=1 !计算Y=F（x） Y(1)=F1(X(1),X(2),X(3),X(4) Y(2)=F2(X(1),X(2),X(3),X(4) Y(3)=F3(X(1),X(2),X(3),X(4) Y(4)=F4(X(1),X(2),X(3),X(4) !计算JA（N,N） E=0.0 DO I=1,N DO J=1,N DO JJ=1,N IF(JJ=J) THEN E(JJ)=X(JJ)+H(JJ) ELSE E(JJ)=X(JJ) ENDIF ENDDO IF(I=1) THEN JA(I,J)=(F1(E(1),E(2),E(3),E(4)-F1(X(1),X(2),X(3),X(4)/H

9、(J) ELSEIF(I=2) THEN JA(I,J)=(F2(E(1),E(2),E(3),E(4)-F2(X(1),X(2),X(3),X(4)/H(J) ELSEIF(I=3) THEN JA(I,J)=(F3(E(1),E(2),E(3),E(4)-F3(X(1),X(2),X(3),X(4)/H(J) ELSEIF(I=4) THEN JA(I,J)=(F4(E(1),E(2),E(3),E(4)-F4(X(1),X(2),X(3),X(4)/H(J) ENDIF ENDDO ENDDO !求解线性方程组 CALL GAUSS(JA,XK,-Y,N) !判断精度 CALL NORM

10、(XK,N,E1) CALL NORM(X,N,E2) IF(E1/E2TA).OR.(A(L,K)=TA) THEN TA=A(L,K) TL=L DO J=K,N T(K,J)=A(K,J) A(K,J)=A(TL,J) A(TL,J)=T(K,J) ENDDO TB(K)=B(K) B(K)=B(TL) B(TL)=TB(K) ENDIF ENDDO DO I=K+1,N M(I,K)=A(I,K)/A(K,K) A(I,K)=0 DO J=K+1,N A(I,J)=A(I,J)-M(I,K)*A(K,J) ENDDO B(I)=B(I)-M(I,K)*B(K) ENDDOENDDO!回

11、代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDO RETURNEND!*求向量的无穷数*SUBROUTINE NORM(X,N,A)DIMENSION X(N)REAL(8) X,AA=ABS(X(1)DO I=2,N IF(ABS(X(I)ABS(X(I-1) THEN A=ABS(X(I) ENDIFENDDORETURNEND!*分片二次代数插值*SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIM

12、ENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!*数据赋值*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,-

13、0.5,-0.66,-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!*计算K,R*IF(UX(N-1)-H/2) THEN K=N-1 ELSE DO I=3,N-2 IF(UX(I)-H/2).AND.(U=X(I)+H/2) THEN K=I ENDIF ENDDOENDIFIF(VY(M-1)-T/2) THEN R=M-1 ELSE DO J=3,M-2 IF(VY(J)-T/2).AND.(VN) P=NIF(P20) P=20IF(QM) Q=MIF(Q20)

14、Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,N APX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,M V(1,J)=0.0 DO I=1,N V(1,J)=V(1,J)+Z(I,J) ENDDO V(1,J)=V(1,J)/D1ENDDOIF(P1) THEN D2=0.0 APX(2)=0.0 DO I=1,N G=X(I)-APX(1) D2=D2+G*G APX(2)=APX(2)+(X(I)-XX)*G*G ENDDO APX(2)=APX(2)/D2 BX(2)=D2/D1 DO J=1,M V(2,J)=0.0 DO

15、I=1,N G=X(I)-APX(1) V(2,J)=V(2,J)+Z(I,J)*G ENDDO V(2,J)=V(2,J)/D2 ENDDO D1=D2ENDIFDO K=3,P D2=0.0 APX(K)=0.0 DO J=1,M V(K,J)=0.0 ENDDO DO I=1,N G1=1.0 G2=X(I)-APX(1) DO J=3,K G=(X(I)-APX(J-1)*G2-BX(J-1)*G1 G1=G2 G2=G ENDDO D2=D2+G*G APX(K)=APX(K)+X(I)*G*G DO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDO ENDDO

16、DO J=1,M V(K,J)=V(K,J)/D2 ENDDO APX(K)=APX(K)/D2 BX(K)=D2/D1 D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,M APY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,P U(J,1)=0.0 DO I=1,M U(J,1)=U(J,1)+V(J,I) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q1)THEN D2=0.0 APY(2)=0.0 DO I=1,M G=Y(I)-APY(1) D2=D2+G*G APY(2)=APY(2)+(Y(I)*G*G ENDD

17、O APY(2)=APY(2)/D2 BY(2)=D2/D1 DO J=1,P U(J,2)=0.0 DO I=1,M G=Y(I)-APY(1) U(J,2)=U(J,2)+V(J,I)*G ENDDO U(J,2)=U(J,2)/D2 ENDDO D1=D2ENDIFDO K=3,Q D2=0.0 APY(K)=0.0 DO J=1,P U(J,K)=0.0 ENDDO DO I=1,M G1=1.0 G2=Y(I)-APY(1) DO J=3,K G=(Y(I)-APY(J-1)*G2-BY(J-1)*G1 G1=G2 G2=G ENDDO D2=D2+G*G APY(K)=APY(K)

18、+Y(I)*G*G DO J=1,P U(J,K)=U(J,K)+V(J,I)*G ENDDO ENDDO DO J=1,P U(J,K)=U(J,K)/D2 ENDDO APY(K)=APY(K)/D2 BY(K)=D2/D1 D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,P DO J=1,Q A(I,J)=0.0 ENDDOENDDODO I=3,Q V(I,I)=V(I-1,I-1) V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2) IF(I=4) THEN DO K=I-2,2,-1 V(I,K)=

19、-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDO ENDIF V(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDODO I=1,P IF(I=1) THEN T(1)=1.0 T1(1)=1.0 ELSEIF(I=2) THEN T(1)=-APX(1) T(2)=1.0 T2(1)=T(1) T2(2)=T(2) ELSE T(I)=T2(I-1) T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2) IF(I=4) THEN DO K=I-2,2,-1 T(K)=-APX(I-1)

20、*T2(K)+T2(K-1)-BX(I-1)*T1(K) ENDDO ENDIF T(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1) T2(I)=T(I) DO K=I-1,1,-1 T1(K)=T2(K) T2(K)=T(K) ENDDO ENDIF DO J=1,Q DO K=I,1,-1 DO L=J,1,-1 A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L) ENDDO ENDDO ENDDOENDDODT1=0.0DO I=1,N X1=X(I) DO J=1,M Y1=Y(J) X2=1.0 DD=0.0 DO K=1,P G=A(K,Q) DO

21、KK=Q-1,1,-1 G=G*Y1+A(K,KK) ENDDO G=G*X2 DD=DD+G X2=X2*X1 ENDDO DT=DD-Z(I,J) DT1=DT1+DT*DT ENDDOENDDO RETURNEND 三、计算结果数表(x,y,f(x,y): X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+00 0.00 0.550 1.322 0.269 0.66E+00 0.00 0.600 1.299 0.295 0.35E+00 0.00 0.650 1.277 0.322 0.94E+00 0.00 0.700 1.255 0.350

22、0.30E-02 0.00 0.750 1.235 0.377 -0.87E-01 0.00 0.800 1.215 0.406 -0.58E+00 0.00 0.850 1.196 0.434 -0.72E+00 0.00 0.900 1.177 0.463 -0.54E+00 0.00 0.950 1.159 0.492 -0.86E+00 0.00 1.000 1.142 0.521 -0.01E+00 0.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.0

23、0E+00 0.00 1.200 1.079 0.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.300 1.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.00 1.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+00 0.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00

24、 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+00 0.08 1.350 1.106 0.709 -0.98E+00 0.08 1.400 1.0