diff --git a/DG-BKG/ARK/src/DG2D_kinetic.f90~ b/DG-BKG/ARK/src/DG2D_kinetic.f90~ deleted file mode 100644 index 8733fc6..0000000 --- a/DG-BKG/ARK/src/DG2D_kinetic.f90~ +++ /dev/null @@ -1,154 +0,0 @@ -Program DG2D_kinetic -use universal_const ! universal.f90 -use Legendre ! Legendre.f90 -use RK_Var -use MD2D_Grid -use State_Var -use Metric_Var -use NorVec_Var -use Kinetic_Var -use Char_Var -implicit none - -!real(kind=8) :: time_start,time_end,sqTDM,ttau_scale,dx,CCFL -!real(kind=8) :: a_vec(1:2), -real(kind=8) :: error_L2,error_max,dtt,dt_eq -real(kind=8) :: a_max -integer :: ierr,q,lid,rks, sqrTDM, sqrTDM2,sqrTDM4, sqrTDM8 - -! --- Program Begin --- -! --- Initialize stage --- -! --- Initialize universal constants --- - - call init_universal_const !-- universal.f90 - - call Input_Domain_Grid_Parameters !-- MD2D_Grid.f90 - - RK_Stage=6 - - lid=22 - open(22,file='./Parameter/in.in',form='formatted', status='old') - read(lid,*) !==============================================! - read(lid,*) Time_final - read(lid,*) a_vec(1) - read(lid,*) a_vec(2) - read(lid,*) !==============================================! - -! --- Initialize Domain Parameter --- -! --- Initialize Legendre Pseudospectral Module --- -! --- 1. Gauss-Lobatto-Legendre Grid --- -! --- 2. Gauss-Lobatto-Legendre Quadurature Weights --- -! --- 3. Legenedre Differentiation Matrix --- - - call alloc_mem_u(PDeg1,PDeg2,PND1,PND2,TotNum_DM,RK_Stage) !-- Legendre.f90 - - call Init_LGL_Diff_Matrix_23D(PND1,PDeg1) !-- Legendre.f90 - - call Init_Physical_Grid_2D !-- Grid2D_Pack.f90 - - call Init_intept_Grid_2D !-- Grid2D_Pack.f90 - - call Init_Metric !-- Metric_Pack.f90 - - call Init_Normal_Vector !-- NorVec_Pack.f90 - - call alloc_mem_Char_Var !-- Char_Var.f90 - - call Init_Penalty_Parameters(ttau_scale) !-- BC_Pack.f90 - - call Init_Kinetic_Variables(PolyNodes,TotNum_DM) !-- Kinetic_Var - -! --- set up initial condition in DG2D_Initial --- -! --- Initialize Time Integration --- - -! call DG2D_Initial(a_vec,Time_final) !-- NorVec_Pack.f90 -! call DG2D_Initial(Time_final) !-- NorVec_Pack.f90 - call DG2D_Kinetic_Initial(Time_final,dt_eq) !-- Kinetic_Pack.f90 - - a_max=maxval(GH) - write(6,*) "a_max",a_max - sqTDM = sqrt(dble(TotNum_DM)) - - sqrTDM=nint(sqTDM) - sqrTDM2=sqrTDM/2 - sqrTDM4=sqrTDM/4 - sqrTDM8=sqrTDM/8 - - ttau_scale = dble(2/sqTDM) - dx = 2/sqTDM - CCFL = dble(2*(2*PDeg1+1))*a_max - CCFL = dble(1/CCFL) *cfl_in -!0.1d0 - dtt = dx*CCFL - rks = 6 !PDeg1 + 1 - - -! call DG2D_Edge(a_vec) !-- NorVec_Pack.f90 - call DG2D_Edge() !-- NorVec_Pack.f90 - -! call Init_SSPRK(rk) !-- RK_Var.f90 -! call Init_RK4S5() !-- RK_Var.f90 - call Init_ARK() !-- RK_Var.f90 - - write(*,*)"dt = ",dtt,"Deg=",PDeg1,"dt_eq",dt_eq - -! --- Time advancing --- - - time = 0.0 - call cpu_time(time_start) - - do while (time .lt. Time_final) - if ( (time+dtt) .ge. Time_final) then - dtt = Time_final-time - endif - Fold=F_alt - F_new = F_alt -! F_alt = F_new - - CALL CALDOM - - do q = 1,RK_Stage - - CALL EQUILIBRIUM - - call DG2D_flux_ark(q) ! Diff.Pack.f90 - -! call SSPRK_marching_2D(q,rk,dtt) ! RK_Var.f90 -! call RK4s5_marching_2D_kinetic(q,rk,dtt) ! RK_Var.f90 - call ARK_marching_2D(q,dtt) - -! call Comp_RU(nnx,pp) -! call Comp_ZTP(nnx,pp) - CALL CALDOM !compute variables from F_alt - CALL EQUILIBRIUM - enddo - F_alt = F_new -! F_new = F_alt - -! write(6,*) "Dm 28",x1(PND1,PND2,28),x2(PND1,PND2,28),"Density",R_loc(PND1,PND2,28) -! write(6,*) "Dm 29",x1(0,PND2,29),x2(0,PND2,29),"Density",R_loc(0,PND2,29) -! write(6,*) "Dm 36",x1(PND1,0,36),x2(PND1,0,36),"Density",R_loc(PND1,0,36) -! write(6,*) "Dm 37",x1(0,0,37),x2(0,0,37),"Density",R_loc(0,0,37) - -DDK=sqrTDM*(sqrTDM2-1)+sqrTDM2 -!DDK=sqrTDM*(3*sqrTDM8-1)+sqrTDM8*3 -write(*,777) DDK,x1(PND1,PND2,DDK),x2(PND1,PND2,DDK),R_loc(PND1,PND2,DDK) -write(*,777) DDK+1,x1(0,PND2,DDK+1),x2(0,PND2,DDK+1),R_loc(0,PND2,DDK+1) -write(*,777) DDK+sqrTDM,x1(PND1,0,DDK+sqrTDM),x2(PND1,0,DDK+sqrTDM),R_loc(PND1,0,DDK+sqrTDM) -write(*,777) DDK+sqrTDM+1,x1(0,0,DDK+sqrTDM+1),x2(0,0,DDK+sqrTDM+1),R_loc(0,0,DDK+sqrTDM+1) -777 format (1X,'Domain #:',I7,2X,'at x=:',F7.4,2X,'at y=:',F7.4,2X,'Density=',F7.4) - - time = time + dtt -write(6,*) "time=",time - - enddo ! do while - - call cpu_time(time_end) - - call DG_compute_error(error_L2,error_max,Time_final) !,a_vec) ! Diff.Pack.f90 - call Tecplt_Output(dx*CCFL) - - write(*,*)"error_L2=",error_L2,"error_max=",error_max - write(*,"(a,1x,1f18.6,1x,a)") "Fortran, cost",time_end-time_start,"s" - -end Program DG2D_kinetic \ No newline at end of file diff --git a/DG-BKG/ARK/src/RK_Pack.f90~ b/DG-BKG/ARK/src/RK_Pack.f90~ deleted file mode 100644 index 73b4d60..0000000 --- a/DG-BKG/ARK/src/RK_Pack.f90~ +++ /dev/null @@ -1,108 +0,0 @@ - - !---------------------------------------------------------------------------- - subroutine SSPRK_marching_2D(q,rk,dtt) - use MD2D_Grid - use Legendre -use RK_Var - - implicit none - integer:: rk,q - real(kind=8):: dtt - - if (q .lt. rk) then - F_new = F_new + alpha(q+1)*(F_alt+dtt*F_t) - F_alt = F_alt + dtt*F_t - else - F_new = F_new + alpha(rk)*dtt*F_t - endif -! write(6,*) "F_t",F_t(1:3,1,0:2,0:2,3) - end subroutine SSPRK_marching_2D - !------------------------------------------------------------------------------ - subroutine RK4s5_marching_2D_kinetic(q,rk,dtt) - !--Declare subroutine arguments - use MD2D_Grid - use Legendre -use RK_Var - - implicit none - integer:: rk,q - real(kind=8):: dtt - - - F_new = RK_para_a(q)*F_new+ dtt*F_t - F_alt = F_alt + RK_para_b(q)*F_new - -! q_tmp(0:N1) = RK_para_a(ks)*q_tmp(0:N1) & -! + dt*dqdt(0:N1,0:N2) -! q(0:N1,0:N2) = q(0:N1,0:N2) & -! + RK_para_b(ks)*q_tmp(0:N1,0:N2) - - - return - ! - end subroutine RK4s5_marching_2D_kinetic - -!===================================================================== - - subroutine ERK_marching_2D(l,dt_rk) - use State_Var - use RK_Var - use Legendre - implicit none -! integer:: nnx ,pp ,l,K - integer:: l - real(kind=8):: dt_rk -! Local Variables - integer:: i,j - -! F_new F at the next time-step <=> F_new(1D) -! F_alt F at the new stage F_alt(2D) <=> F(1D) - - if (l .LT. RK_Stage) then - - F_new=F_new+ dt_rk*const_b(l)*(F_s(:,:,:,:,:,l)+F_ns(:,:,:,:,:,l)) - F_alt=Fold - do j=1,l - F_alt = F_alt + dt_rk*(const_a_I(l+1,j)*F_s(:,:,:,:,:,j) + & - const_a_E(l+1,j)*F_ns(:,:,:,:,:,j)) - enddo - else - F_new=F_new+ dt_rk*const_b(l)*(F_s(:,:,:,:,:,l)+F_ns(:,:,:,:,:,l)) - endif - -! call Comp_RU(nnx,pp) -! call Comp_ZTP(nnx,pp) - end subroutine ERK_marching_2D - - !------------------------------------------------------------------------------ - subroutine ARK_marching_2D(l,dt_rk) - use State_Var - use RK_Var - use Legendre - implicit none -! integer:: nnx ,pp ,l,K - integer:: l - real(kind=8):: dt_rk -! Local Variables - integer:: i,j - -! F_new F at the next time-step <=> F_new(1D) -! F_alt F at the new stage F_alt(2D) <=> F(1D) - - if (l .LT. RK_Stage) then - - F_new=F_new+ dt_rk*const_b(l)*(F_s(:,:,:,:,:,l)+F_ns(:,:,:,:,:,l)) - F_alt=Fold - do j=1,l - F_alt = F_alt + dt_rk*(const_a_I(l+1,j)*F_s(:,:,:,:,:,j) + & - const_a_E(l+1,j)*F_ns(:,:,:,:,:,j)) - enddo - else - F_new=F_new+ dt_rk*const_b(l)*(F_s(:,:,:,:,:,l)+F_ns(:,:,:,:,:,l)) - endif - -! call Comp_RU(nnx,pp) -! call Comp_ZTP(nnx,pp) - end subroutine ARK_marching_2D - - !------------------------------------------------------------------------------ diff --git a/myf90/Daniel_Furano/LongManual.pdf b/myf90/Daniel_Furano/LongManual.pdf new file mode 100644 index 0000000..c3454fe Binary files /dev/null and b/myf90/Daniel_Furano/LongManual.pdf differ diff --git a/myf90/Daniel_Furano/ShortManual.pdf b/myf90/Daniel_Furano/ShortManual.pdf new file mode 100644 index 0000000..f9fef4d Binary files /dev/null and b/myf90/Daniel_Furano/ShortManual.pdf differ diff --git a/myf90/Daniel_Furano/splib.txt b/myf90/Daniel_Furano/splib.txt new file mode 100644 index 0000000..5b2b8ee --- /dev/null +++ b/myf90/Daniel_Furano/splib.txt @@ -0,0 +1,3549 @@ +***************************************************************** +* FORTRAN ROUTINES FOR SPECTRAL METHODS +* BY DANIELE FUNARO +* DEPARTMENT OF MATHEMATICS +* UNIVERSITY OF MODENA +* VIA CAMPI 213/B, 41100 MODENA, ITALY +* E-MAIL: funaro@unimo.it +***************************************************************** +C + SUBROUTINE GAMMAF(X,GX) +***************************************************************** +* COMPUTES THE GAMMA FUNCTION AT A GIVEN POINT +* X = ARGUMENT GREATER THAN 1.E-75 AND SMALLER THAN 57. +* GX= VALUE OF GAMMA IN X +***************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION C(11) + + IF (X .LE. 1.D-75 .OR. X .GE. 57.D0) THEN + WRITE(*,*) 'ARGUMENT OUT OF RANGE IN SUBROUTINE GAMMAF' + RETURN + ENDIF + + PI = 3.14159265358979323846D0 + EPS = 1.D-14 + XX = X + GX = 1.0D0 + +1 IF (DABS(XX-1.D0) .LT. EPS) RETURN + IF (XX .GE. 1.D0) THEN + XX = XX-1.D0 + GX = GX*XX + GOTO 1 + ENDIF + IND = 0 + IF (XX .LT. .5D0) THEN + IND = 1 + GX = GX*PI/DSIN(PI*XX) + XX = 1.D0-XX + ENDIF + + PR = 1.D0 + S = 0.426401432711220868D0 + C(1) = -0.524741987629368444D0 + C(2) = 0.116154405493589130D0 + C(3) = -0.765978624506602380D-2 + C(4) = 0.899719449391378898D-4 + C(5) = -0.194536980009534621D-7 + C(6) = 0.199382839513630987D-10 + C(7) = -0.204209590209541319D-11 + C(8) = 0.863896817907000175D-13 + C(9) = 0.152237501608472336D-13 + C(10) = -0.82572517527771995D-14 + C(11) = 0.29973478220522461D-14 + + DO 2 K=1,11 + PR = PR*(XX-DFLOAT(K))/(XX+DFLOAT(K-1)) + S = S+C(K)*PR +2 CONTINUE + G = S*DEXP(1.D0-XX)*(XX+4.5D0)**(XX-.5D0) + IF (IND .EQ. 1) THEN + GX = GX/G + ELSE + GX = GX*G + ENDIF + + RETURN + END +C + SUBROUTINE VAJAPO(N,A,B,X,Y,DY,D2Y) +************************************************************ +* COMPUTES THE VALUE OF THE JACOBI POLYNOMIAL OF DEGREE N +* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT +* N = DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + AB = A+B + Y = .5D0*(AB+2.D0)*X+.5D0*(A-B) + DY = .5D0*(AB+2.D0) + D2Y = 0.D0 + IF(N .EQ. 1) RETURN + + YP = 1.D0 + DYP = 0.D0 + D2YP = 0.D0 + DO 1 I=2,N + DI = DFLOAT(I) + C0 = 2.D0*DI+AB + C1 = 2.D0*DI*(DI+AB)*(C0-2.D0) + C2 = (C0-1.D0)*(C0-2.D0)*C0 + C3 = (C0-1.D0)*(A-B)*AB + C4 = 2.D0*(DI+A-1.D0)*C0*(DI+B-1.D0) + YM = Y + Y = ((C2*X+C3)*Y-C4*YP)/C1 + YP = YM + DYM = DY + DY = ((C2*X+C3)*DY-C4*DYP+C2*YP)/C1 + DYP = DYM + D2YM = D2Y + D2Y = ((C2*X+C3)*D2Y-C4*D2YP+2.D0*C2*DYP)/C1 + D2YP = D2YM +1 CONTINUE + + RETURN + END +C + SUBROUTINE VALEPO(N,X,Y,DY,D2Y) +************************************************************** +* COMPUTES THE VALUE OF THE LEGENDRE POLYNOMIAL OF DEGREE N +* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT +* N = DEGREE OF THE POLYNOMIAL +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + Y = X + DY = 1.D0 + D2Y = 0.D0 + IF(N .EQ. 1) RETURN + + YP = 1.D0 + DYP = 0.D0 + D2YP = 0.D0 + DO 1 I=2,N + C1 = DFLOAT(I) + C2 = 2.D0*C1-1.D0 + C4 = C1-1.D0 + YM = Y + Y = (C2*X*Y-C4*YP)/C1 + YP = YM + DYM = DY + DY = (C2*X*DY-C4*DYP+C2*YP)/C1 + DYP = DYM + D2YM = D2Y + D2Y = (C2*X*D2Y-C4*D2YP+2.D0*C2*DYP)/C1 + D2YP = D2YM +1 CONTINUE + + RETURN + END +C + SUBROUTINE VACHPO(N,X,Y,DY,D2Y) +*************************************************************** +* COMPUTES THE VALUE OF THE CHEBYSHEV POLYNOMIAL OF DEGREE N +* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT +* N = DEGREE OF THE POLYNOMIAL +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +*************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + Y = X + DY = 1.D0 + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + YP = 1.D0 + DYP = 0.D0 + D2YP = 0.D0 + DO 1 K=2,N + YM = Y + Y = 2.D0*X*Y-YP + YP = YM + DYM = DY + DY = 2.D0*X*DY+2.D0*YP-DYP + DYP = DYM + D2YM= D2Y + D2Y = 2.D0*X*D2Y+4.D0*DYP-D2YP + D2YP= D2YM +1 CONTINUE + + RETURN + END +C + SUBROUTINE VALAPO(N,A,X,Y,DY,D2Y) +************************************************************** +* COMPUTES THE VALUE OF THE LAGUERRE POLYNOMIAL OF DEGREE N +* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT +* N = DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + Y = 1.D0+A-X + DY = -1.D0 + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + YP = 1.D0 + DYP = 0.D0 + D2YP= 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + B1 = (2.D0*DK+A-1.D0-X)/DK + B2 = (DK+A-1.D0)/DK + YM = Y + Y = B1*Y-B2*YP + YP = YM + DYM = DY + DY = B1*DY-YP/DK-B2*DYP + DYP = DYM + D2YM= D2Y + D2Y = B1*D2Y-2.D0*DYP/DK-B2*D2YP + D2YP= D2YM + 1 CONTINUE + + RETURN + END +C + SUBROUTINE VAHEPO(N,X,Y,DY,D2Y) +************************************************************* +* COMPUTES THE VALUE OF THE HERMITE POLYNOMIAL OF DEGREE N +* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT +* N = DEGREE OF THE POLYNOMIAL +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + Y = 2.D0*X + DY = 2.D0 + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + YP=1.D0 + DO 1 K=2,N + DK = DFLOAT(K-1) + YM = Y + Y = 2.D0*X*Y-2.D0*DK*YP + YPM = YP + YP = YM +1 CONTINUE + DN = 2.D0*DFLOAT(N) + DNN = 2.D0*DFLOAT(N-1) + DY = DN*YP + D2Y = DN*DNN*YPM + + RETURN + END +C + SUBROUTINE VALASF(N,A,X,Y,DY) +********************************************************************* +* COMPUTES THE VALUES OF THE SCALED LAGUERRE FUNCTION OF DEGREE N +* AND ITS FIRST DERIVATIVE AT A GIVEN POINT +* N = DEGREE +* A = PARAMETER >-1 +* X = POINT (NON NEGATIVE) IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE FUNCTION IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + IF (N .EQ. 0) RETURN + + C0 = 1.D0/(4.D0+X) + C1 = 4.D0*C0/(A+1.D0) + Y = C1*(A+1.D0-X) + DY = -C0*C1*(A+5.D0) + IF (N .EQ. 1) RETURN + + YP = 1.D0 + DYP = 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + DK4 = 4.D0*DK + C0 = 1.D0/(DK4+X) + C1 = DK4+X-2.D0 + C2 = 1.D0/(C1-2.D0) + C3 = DK4*C0/(DK+A) + C4 = 2.D0*DK+A-1.D0 + C5 = C4-X + C6 = C4+DK4 + C7 = C2*DFLOAT(4*(K-1)**2) + DYM = DY + DY = C3*(C5*DY-C0*C6*Y+C7*(2.D0*C0*C1*C2*YP-DYP)) + DYP = DYM + YM = Y + Y = C3*(C5*Y-C7*YP) + YP = YM +1 CONTINUE + + RETURN + END +C + SUBROUTINE VAHESF(N,X,Y,DY) +*********************************************************** +* COMPUTES THE VALUES OF THE SCALED HERMITE FUNCTION OF +* DEGREE N AND ITS FIRST DERIVATIVE AT A GIVEN POINT +* N = DEGREE +* X = POINT IN WHICH THE COMPUTATION IS PERFORMED +* Y = VALUE OF THE FUNCTION IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +*********************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + Y = 1.D0 + DY = 0.D0 + IF (N .EQ. 0) RETURN + + XX = X*X + M = N/2 + IF(N .EQ. 2*M) THEN + A = -.5D0 + CALL VALASF(M,A,XX,Y,DY) + DY = 2.D0*X*DY + ELSE + A = .5D0 + CALL VALASF(M,A,XX,Y,DY) + DY = Y+2.D0*XX*DY + Y = X*Y + ENDIF + + RETURN + END +C + SUBROUTINE ZEJAGA(N,A,B,CS,DZ) +*************************************************************** +* COMPUTES THE ZEROES OF THE JACOBI POLYNOMIAL OF DEGREE N +* N = THE NUMBER OF ZEROES +* A = PARAMETER BETWEEN -1/2 AND 1/2 +* B = PARAMETER BETWEEN -1/2 AND 1/2 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +*************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1) + IF (N .EQ .0) RETURN + + AB = A+B + CS(1) = (B-A)/(AB+2.D0) + DZ(1) = .5D0*AB+1.D0 + IF(N .EQ. 1) RETURN + + EPS= 1.E-17 + PH = 1.57079632679489661923D0 + C = PH/(2.D0*DFLOAT(N)+AB+1.D0) + DO 1 I=1,N + DI = DFLOAT(I) + CSX = -DCOS(C*(4.D0*DI+AB-1.D0)) + DO 2 IT=1,8 + CALL VAJAPO(N,A,B,CSX,Y,DY,D2Y) + IF(DABS(Y) .LT. EPS) GOTO 3 + CSX = CSX-Y/DY +2 CONTINUE +3 IF(DABS(CSX) .LT. EPS) CSX=0.D0 + CS(I) = CSX + DZ(I) = DY +1 CONTINUE + + RETURN + END +C + SUBROUTINE ZELEGA(N,CS,DZ) +*************************************************************** +* COMPUTES THE ZEROES OF THE LEGENDRE POLYNOMIAL OF DEGREE N +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +*************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1) + IF (N .EQ .0) RETURN + + CS(1) = 0.D0 + DZ(1) = 1.D0 + IF(N .EQ. 1) RETURN + + N2 = N/2 + IN = 2*N-4*N2-1 + PH = 1.57079632679489661923D0 + C = PH/(2.D0*DFLOAT(N)+1.D0) + DO 1 I=1,N2 + DI = DFLOAT(I) + CSX = DCOS(C*(4.D0*DI-1.D0)) + DO 2 IT=1,8 + CALL VALEPO(N,CSX,Y,DY,D2Y) + CSX = CSX-Y/DY +2 CONTINUE + CS(I) = -CSX + CS(N-I+1) = CSX + DZ(I) = DY*DFLOAT(IN) + DZ(N-I+1) = DY +1 CONTINUE + + IF(IN .EQ. -1) RETURN + CSX = 0.D0 + CS(N2+1) = CSX + CALL VALEPO(N,CSX,Y,DY,D2Y) + DZ(N2+1) = DY + + RETURN + END +C + SUBROUTINE ZECHGA(N,CS,DZ) +**************************************************************** +* COMPUTES THE ZEROES OF THE CHEBYSHEV POLYNOMIAL OF DEGREE N +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1) + IF (N .EQ. 0) RETURN + + CS(1) = 0.D0 + DZ(1) = 1.D0 + IF(N .EQ. 1) RETURN + + N2 = N/2 + IN = 1+4*N2-2*N + PH = 1.57079632679489661923D0 + DN = DFLOAT(N) + C = PH/DN + SI = -1.D0 + DO 1 I=1,N2 + DI = DFLOAT(I) + CSX = DCOS(C*(2.D0*DI-1.D0)) + CS(I) = -CSX + CS(N-I+1) = CSX + QX = DN/DSQRT(1.D0-CSX*CSX) + DZ(I) = QX*SI*DFLOAT(IN) + DZ(N-I+1) = -QX*SI + SI = -SI +1 CONTINUE + + IF(IN .EQ. 1) RETURN + CS(N2+1) = 0.D0 + N4 = N2/2 + IN2 = 1+4*N4-2*N2 + DZ(N2+1) = DN*DFLOAT(IN2) + + RETURN + END +C + SUBROUTINE ZELAGA(N,A,CS,DZ) +************************************************************************ +* COMPUTES THE ZEROES OF THE LAGUERRE POLYNOMIAL OF DEGREE N +* N = THE NUMBER OF ZEROES +* A = PARAMETER >-1 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = DERIVATIVES OF THE SCALED FUNCTIONS AT THE ZEROES, DZ(I), I=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + CS(1) = A1 + DZ(1) = -4.D0/(A1*(A+5.D0)) + IF (N .EQ. 1) RETURN + + PI = 3.14159265358979323846D0 + DN = DFLOAT(N) + C1 = 2.D0*DN+A1 + DO 1 M=1,N + DM = DFLOAT(M) + C2 = 2.D0*(DN+.75D0-DM)*PI/C1 + XN = (C2+PI)/2.D0 + DO 2 IT=1,8 + XP = XN + XN = (DSIN(XP)-XP*DCOS(XP)+C2)/(1.D0-DCOS(XP)) +2 CONTINUE + Z = (DCOS(XN/2.D0))**2 + ZD = 1.D0/(Z-1.D0) + CSX= 2.D0*C1*Z-((1.25D0*ZD+1.D0)*ZD-1.D0+3.D0*A*A)/(6.D0*C1) + DO 3 IT=1,6 + CALL VALASF(N,A,CSX,Y,DY) + CSX = CSX-Y/DY +3 CONTINUE + CS(M) = CSX + DZ(M) = DY +1 CONTINUE + +4 IN = 0 + DO 5 M=1,N-1 + IF(CS(M) .LE. CS(M+1)) GOTO 5 + CSM = CS(M) + CS(M) = CS(M+1) + CS(M+1) = CSM + DZM = DZ(M) + DZ(M) = DZ(M+1) + DZ(M+1) = DZM + IN = 1 +5 CONTINUE + IF(IN .EQ. 1) GOTO 4 + + RETURN + END +C + SUBROUTINE ZEHEGA(N,CS,DZ) +************************************************************************ +* COMPUTES THE ZEROES OF THE HERMITE POLYNOMIAL OF DEGREE N +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = DERIVATIVES OF THE SCALED FUNCTIONS AT THE ZEROES, DZ(I), I=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1) + IF (N .EQ. 0) RETURN + + M = N/2 + CS(M+1) = 0.D0 + DZ(M+1) = 1.D0 + IF(N .EQ. 1) RETURN + + IN = 2*N-4*M-1 + IF(IN .EQ. -1) THEN + A = -.5D0 + CALL ZELAGA(M,A,CS,DZ) + ELSE + A = .5D0 + CALL ZELAGA(M,A,CS,DZ) + ENDIF + + DO 1 I=1,M + CSX = DSQRT(CS(M-I+1)) + CS(N-I+1) = CSX + IF(IN .EQ. -1) THEN + DZ(N-I+1) = 2.D0*CSX*DZ(M-I+1) + ELSE + DZ(N-I+1) = 2.D0*CSX*CSX*DZ(M-I+1) + ENDIF +1 CONTINUE + + DO 2 I=1,M + CS(I) = -CS(N-I+1) + DZ(I) = DZ(N-I+1)*DFLOAT(IN) +2 CONTINUE + + RETURN + END +C + SUBROUTINE ZEJAGL(N,A,B,ET,VN) +******************************************************************** +* COMPUTES THE NODES RELATIVE TO THE JACOBI GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER BETWEEN -1/2 AND 1/2 +* B = PARAMETER BETWEEN -1/2 AND 1/2 +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +******************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*) + IF (N .EQ. 0) RETURN + + ET(0) = -1.D0 + ET(N) = 1.D0 + X = -1.D0 + CALL VAJAPO(N,A,B,X,Y,DY,D2Y) + VN(0) = Y + X = 1.D0 + CALL VAJAPO(N,A,B,X,Y,DY,D2Y) + VN(N) = Y + IF (N .EQ. 1) RETURN + + EPS= 1.E-17 + PH = 1.57079632679489661923D0 + AB = A+B + DN = DFLOAT(N) + C = PH/(2.D0*DN+AB+1.D0) + N1 = N-1 + A1 = A+1.D0 + B1 = B+1.D0 + DO 1 I=1,N1 + DI = DFLOAT(I) + ETX = -DCOS(C*(4.D0*DI+AB+1.D0)) + DO 2 IT=1,8 + CALL VAJAPO(N1,A1,B1,ETX,Y,DY,D2Y) + IF(DABS(Y) .LE. EPS) GOTO 3 + ETX = ETX-Y/DY +2 CONTINUE +3 IF(DABS(ETX) .LE. EPS) ETX=0.D0 + ET(I) = ETX + VN(I) = -.5D0*DY*(1.D0-ETX*ETX)/DN +1 CONTINUE + + RETURN + END +C + SUBROUTINE ZELEGL(N,ET,VN) +********************************************************************* +* COMPUTES THE NODES RELATIVE TO THE LEGENDRE GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*) + IF (N .EQ. 0) RETURN + + N2 = (N-1)/2 + SN = DFLOAT(2*N-4*N2-3) + ET(0) = -1.D0 + ET(N) = 1.D0 + VN(0) = SN + VN(N) = 1.D0 + IF (N .EQ. 1) RETURN + + ET(N2+1) = 0.D0 + X = 0.D0 + CALL VALEPO(N,X,Y,DY,D2Y) + VN(N2+1) = Y + IF(N .EQ. 2) RETURN + + PI = 3.14159265358979323846D0 + C = PI/DFLOAT(N) + DO 1 I=1,N2 + ETX = DCOS(C*DFLOAT(I)) + DO 2 IT=1,8 + CALL VALEPO(N,ETX,Y,DY,D2Y) + ETX = ETX-DY/D2Y +2 CONTINUE + ET(I) = -ETX + ET(N-I) = ETX + VN(I) = Y*SN + VN(N-I) = Y +1 CONTINUE + + RETURN + END +C + SUBROUTINE ZECHGL(N,ET) +********************************************************************** +* COMPUTES THE NODES RELATIVE TO THE CHEBYSHEV GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* ET = VECTOR OF THE NODES, ET(I), I=0,N +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*) + IF (N .EQ. 0) RETURN + + ET(0) = -1.D0 + ET(N) = 1.D0 + IF (N .EQ. 1) RETURN + + N2 = (N-1)/2 + ET(N2+1) = 0.D0 + IF(N .EQ. 2) RETURN + + PI = 3.14159265358979323846D0 + C = PI/DFLOAT(N) + DO 1 I=1,N2 + ETX = DCOS(C*DFLOAT(I)) + ET(I) = -ETX + ET(N-I) = ETX +1 CONTINUE + + RETURN + END +C + SUBROUTINE ZELAGR(N,A,ET,VN) +************************************************************************ +* COMPUTES THE NODES RELATIVE TO THE LAGUERRE GAUSS-RADAU FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), 0=1,N-1 +* VN = SCALED LAGUERRE FUNCTION AT THE NODES, VN(I), I=0,N-1 +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*) + IF (N .EQ. 0) RETURN + + ET(0) = 0.D0 + VN(0) = 1.D0 + IF (N .EQ. 1) RETURN + + A1 = A+1.D0 + PI = 3.14159265358979323846D0 + N1 = N-1 + DN = DFLOAT(N1) + C1 = 2.D0*DN+A1+1.D0 + DO 1 M=1,N1 + DM = DFLOAT(M) + C2 = 2.D0*(DN+.75D0-DM)*PI/C1 + XN = (C2+PI)/2.D0 + DO 2 IT=1,8 + XP = XN + XN = (DSIN(XP)-XP*DCOS(XP)+C2)/(1.D0-DCOS(XP)) +2 CONTINUE + Z = (DCOS(XN/2.D0))**2 + ZD = 1.D0/(Z-1.D0) + ETX= 2.D0*C1*Z-((1.25D0*ZD+1.D0)*ZD-1.D0+3.D0*A1*A1)/(6.D0*C1) + DO 3 IT=1,6 + CALL VALASF(N1,A1,ETX,Y,DY) + ETX = ETX-Y/DY +3 CONTINUE + ET(M) = ETX + CALL VALASF(N,A,ETX,Y,DY) + VN(M) = Y +1 CONTINUE + +4 IN = 0 + DO 5 M=1,N-2 + IF(ET(M) .LE. ET(M+1)) GOTO 5 + ETM = ET(M) + ET(M) = ET(M+1) + ET(M+1) = ETM + VNM = VN(M) + VN(M) = VN(M+1) + VN(M+1) = VNM + IN = 1 +5 CONTINUE + IF(IN .EQ. 1) GOTO 4 + + RETURN + END +C + SUBROUTINE WEJAGA(N,A,B,CS,DZ,WE) +**************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE JACOBI GAUSS FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER > -1 +* B = PARAMETER > -1 +* CS = ZEROES OF THE JACOBI POLYNOMIAL, CS(I), I=1,N +* DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), WE(1) + IF (N .EQ. 0) RETURN + + AB = A+B+2.D0 + A2 = A+2.D0 + B2 = B+2.D0 + CALL GAMMAF(A2,GA2) + CALL GAMMAF(B2,GB2) + CALL GAMMAF(AB,GAB) + C = .5D0*(2.D0**AB)*GA2*GB2/GAB + DO 1 M=2,N + DM = DFLOAT(M) + C = C*(DM+A)*(DM+B)/(DM*(DM+A+B)) +1 CONTINUE + + DO 2 I=1,N + X = CS(I) + DY = DZ(I) + WE(I) = C/((1.D0-X*X)*DY*DY) +2 CONTINUE + + RETURN + END +C + SUBROUTINE WELEGA(N,CS,DZ,WE) +***************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE LEGENDRE GAUSS FORMULA +* N = ORDER OF THE FORMULA +* CS = ZEROES OF THE LEGENDRE POLYNOMIAL, CS(I), I=1,N +* DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +***************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), WE(1) + IF (N .EQ. 0) RETURN + + N2 = N/2 + DO 1 I=1,N2 + X = CS(I) + DY = DZ(I) + WEX = 2.D0/((1.D0-X*X)*DY*DY) + WE(I) = WEX + WE(N-I+1) = WEX +1 CONTINUE + + IF(N .EQ. 2*N2) RETURN + DY = DZ(N2+1) + WE(N2+1) = 2.D0/(DY*DY) + + RETURN + END +C + SUBROUTINE WECHGA(N,WE) +***************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE CHEBYSHEV GAUSS FORMULA +* N = ORDER OF THE FORMULA +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +***************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION WE(1) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + C = PI/DFLOAT(N) + DO 1 I=1,N + WE(I) = C +1 CONTINUE + + RETURN + END +C + SUBROUTINE WELAGA(N,A,CS,WE) +**************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE LAGUERRE GAUSS FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER > -1 +* CS = ZEROES OF THE LAGUERRE POLYNOMIAL, CS(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), WE(1) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + CALL GAMMAF(A1,GA1) + N1 = N+1 + DN = DFLOAT(N1) + C = GA1/DN + DO 1 M=1,N + DM = DFLOAT(M) + C = C*(DM+A)/(DM+1.D0) +1 CONTINUE + + DO 2 I=1,N + X =CS(I) + CALL VALAPO(N1,A,X,Y,DY,D2Y) + WE(I) = C*X/(Y*Y) +2 CONTINUE + + RETURN + END +C + SUBROUTINE WEHEGA(N,CS,WE) +**************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE HERMITE GAUSS FORMULA +* N = ORDER OF THE FORMULA +* CS = ZEROES OF THE HERMITE POLYNOMIAL, CS(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), WE(1) + IF (N .EQ. 0) RETURN + + PR = 1.77245385090551588D0 + R2 = 1.41421356237309515D0 + C = PR/DFLOAT(N) + N2 = N/2 + DO 1 I=1,N2 + X = CS(I) + YP = 1.D0 + Y = R2*X + DO 2 K=2,N-1 + DK = DFLOAT(K) + RK = DSQRT(DK) + QK = DSQRT(DK-1.D0) + YM = Y + Y = (R2*X*Y-QK*YP)/RK + YP = YM +2 CONTINUE + WEX = C/(Y*Y) + WE(I) = WEX + WE(N-I+1) = WEX +1 CONTINUE + + IF(N .EQ. 2*N2) RETURN + Y = 1.D0 + DO 3 K=2,N-1,2 + DK = DFLOAT(K) + Y = Y*DSQRT((DK-1.D0)/DK) +3 CONTINUE + WE(N2+1) = C/(Y*Y) + + RETURN + END +C + SUBROUTINE WEJAGL(N,A,B,ET,WT) +********************************************************************* +* COMPUTES THE WEIGHTS RELATIVE TO THE JACOBI GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER > -1 +* B = PARAMETER > -1 +* ET = JACOBI GAUSS-LOBATTO NODES, ET(I), I=0,N +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), WT(0:*) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + AB2 = A1+B1 + CALL GAMMAF(A1,GA1) + CALL GAMMAF(B1,GB1) + CALL GAMMAF(AB2,GAB2) + C = (2.D0**AB)*GA1*GB1/GAB2 + WT(0) = 2.D0*C*A1/AB2 + WT(N) = 2.D0*C*B1/AB2 + IF (N .EQ. 1) RETURN + + N1 = N-1 + DN = DFLOAT(N) + C = C*(2.D0*DN+AB)/(DN+AB+1.D0) + C1 = C*A1/((B+2.D0)*AB2) + C2 = C*B1/((A+2.D0)*AB2) + C3 = .5D0*C*A1*B1 + DO 1 K=1,N-2 + DK = DFLOAT(K) + C1 = C1*(DK+A1)*DK/((DK+AB2)*(DK+B+2.D0)) + C2 = C2*(DK+B1)*DK/((DK+AB2)*(DK+A+2.D0)) + C3 = C3*(DK+A1)*(DK+B1)/((DK+2.D0)*(DK+AB+1.D0)) +1 CONTINUE + + SU = 0.D0 + DO 2 M=1,N1 + SU = SU+ET(M) +2 CONTINUE + WT(0) = C1*(DN-1.D0-SU) + WT(N) = C2*(DN-1.D0+SU) + DO 3 I=1,N1 + X = ET(I) + CALL VAJAPO(N,A,B,X,Y,DY,D2Y) + C4 = -C3/Y + CALL VAJAPO(N1,A,B,X,Y,DY,D2Y) + WT(I) = C4/DY +3 CONTINUE + + RETURN + END +C + SUBROUTINE WELEGL(N,ET,VN,WT) +*********************************************************************** +* COMPUTES THE WEIGHTS RELATIVE TO THE LEGENDRE GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* ET = JACOBI GAUSS-LOBATTO NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N +*********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), WT(0:*) + IF (N .EQ. 0) RETURN + + N2 = (N-1)/2 + DN = DFLOAT(N) + C = 2.D0/(DN*(DN+1.D0)) + DO 1 I=0,N2 + X = ET(I) + Y = VN(I) + WTX = C/(Y*Y) + WT(I) = WTX + WT(N-I) = WTX +1 CONTINUE + + IF(N-1 .EQ. 2*N2) RETURN + X = 0.D0 + Y = VN(N2+1) + WT(N2+1) = C/(Y*Y) + + RETURN + END +C + SUBROUTINE WECHGL(N,WT) +************************************************************************ +* COMPUTES THE WEIGHTS RELATIVE TO THE CHEBYSHEV GAUSS-LOBATTO FORMULA +* N = ORDER OF THE FORMULA +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION WT(0:*) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + C = PI/DFLOAT(N) + C2 = .5D0*C + WT(0) =C2 + WT(N) =C2 + IF (N .EQ. 1) RETURN + + DO 1 I=1,N-1 + WT(I) = C +1 CONTINUE + + RETURN + END +C + SUBROUTINE WELAGR(N,A,ET,WT) +********************************************************************* +* COMPUTES THE WEIGHTS RELATIVE TO THE LAGUERRE GAUSS-RADAU FORMULA +* N = ORDER OF THE FORMULA +* A = PARAMETER > -1 +* ET = LAGUERRE GAUSS-RADAU NODES, ET(I), I=0,N-1 +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N-1 +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), WT(0:*) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + CALL GAMMAF(A1,GA1) + C1 = GA1 + WT(0) = C1 + IF (N .EQ. 1) RETURN + + N1 = N-1 + C2 = GA1 + DO 1 K=1,N1 + DK = DFLOAT(K) + C1 = C1*DK/(DK+A1) + C2 = C2*(DK+A)/(DK+1.D0) +1 CONTINUE + WT(0) = C1 + DO 2 I=1,N1 + X = ET(I) + CALL VALAPO(N,A,X,Y,DY,D2Y) + C3 = C2/Y + CALL VALAPO(N1,A,X,Y,DY,D2Y) + WT(I) = C3/DY +2 CONTINUE + + RETURN + END +C + SUBROUTINE WECHCC(N,WK) +********************************************************************* +* COMPUTES THE WEIGHTS OF THE CLENSHAW-CURTIS FORMULA OF ORDER 2*N +* N = INTEGER PARAMETER +* WK = VECTOR OF THE WEIGHTS, WK(I), I=0,2*N +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION WK(0:*) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + M = 2*N + DN = DFLOAT(N) + DM = DFLOAT(M) + C = 1.D0/(DM*DM-1.D0) + WK(0) = C + WK(M) = C + WK(N) = 1.33333333333333333333D0 + IF (N .EQ. 1) RETURN + + DO 1 J=1,N-1 + DJ = DFLOAT(J) + SU = 1.D0-((-1.D0)**J)*C + DO 2 K=1,N-1 + DK = 2.D0*DFLOAT(K) + SU = SU+2.D0*DCOS(DJ*DK*PI/DM)/(1.D0-DK*DK) +2 CONTINUE + WK(J) = SU/DN + WK(M-J) = SU/DN +1 CONTINUE + + SU = 1.D0-((-1.D0)**N)*C + DO 3 K=1,N-1 + DK = 2.D0*DFLOAT(K) + SU = SU+2.D0*((-1.D0)**K)/(1.D0-DK*DK) +3 CONTINUE + WK(N) = SU/DN + + RETURN + END +C + SUBROUTINE INJAGA(N,A,B,CS,DZ,QZ,X,QX) +************************************************************************ +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE ZEROES OF THE JACOBI POLYNOMIAL +* N = THE NUMBER OF ZEROES +* A = PARAMETER > -1 +* B = PARAMETER > -1 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = JACOBI POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), QZ(1) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + CALL VAJAPO(N,A,B,X,Y,DY,D2Y) + QX = 0.D0 + DO 1 J=1,N + ED = X-CS(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QZ(J) + RETURN + ELSE + QX = QX+QZ(J)*Y/(DZ(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INLEGA(N,CS,DZ,QZ,X,QX) +************************************************************************ +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE ZEROES OF THE LEGENDRE POLYNOMIAL +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = LEGENDRE POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), QZ(1) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + CALL VALEPO(N,X,Y,DY,D2Y) + QX = 0.D0 + DO 1 J=1,N + ED = X-CS(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QZ(J) + RETURN + ELSE + QX = QX+QZ(J)*Y/(DZ(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INCHGA(N,DZ,QZ,X,QX) +************************************************************************ +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE ZEROES OF THE CHEBYSHEV POLYNOMIAL +* N = THE NUMBER OF ZEROES +* DZ = CHEBYSHEV POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION DZ(1), QZ(1) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + PH = 1.57079632679489661923D0 + DN = DFLOAT(N) + C = PH/DN + CALL VACHPO(N,X,Y,DY,D2Y) + QX = 0.D0 + DO 1 J=1,N + ED = X+DCOS(C*(2.D0*DFLOAT(J)-1.D0)) + IF(DABS(ED) .LT. EPS) THEN + QX = QZ(J) + RETURN + ELSE + QX = QX+QZ(J)*Y/(DZ(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INLAGA(N,A,CS,DZ,QZ,X,QX) +********************************************************************* +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE ZEROES OF THE LAGUERRE POLYNOMIAL +* N = THE NUMBER OF ZEROES +* A = PARAMETER > -1 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = SCALED LAGUERRE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), QZ(1) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + CALL VALASF(N,A,X,Y,DY) + QX = 0.D0 + DO 1 J=1,N + ED = X-CS(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QZ(J) + RETURN + ELSE + PR = 1.D0 + DO 2 K=1,N + DK = 4.D0*DFLOAT(K) + PR = PR*(DK+X)/(DK+CS(J)) +2 CONTINUE + QX = QX+QZ(J)*Y*PR/(DZ(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INHEGA(N,CS,DZ,QZ,X,QX) +********************************************************************* +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE ZEROES OF THE HERMITE POLYNOMIAL +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = SCALED HERMITE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), QZ(1) + IF (N .EQ. 0) RETURN + + QX = QZ(1) + IF (N .EQ. 1) RETURN + + EPS = 1.D-14 + CALL VAHESF(N,X,Y,DY) + QX = 0.D0 + DO 1 J=1,N + ED = X-CS(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QZ(J) + RETURN + ELSE + PR = 1.D0 + DO 2 K=1,N/2 + DK = 4.D0*DFLOAT(K) + PR = PR*(DK+X*X)/(DK+CS(J)**2) +2 CONTINUE + QX = QX+QZ(J)*Y*PR/(DZ(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INJAGL(N,A,B,ET,VN,QN,X,QX) +********************************************************************* +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE JACOBI GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER > -1 +* B = PARAMETER > -1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + CALL VAJAPO(N,A,B,X,Y,DY,D2Y) + DN = DFLOAT(N) + C = 1.D0/(DN*(DN+A+B+1.D0)) + QX = QN(0)*C*DY*(B+1.D0)*(X-1.D0)/VN(0) + QX = QX+QN(N)*C*DY*(A+1.D0)*(X+1.D0)/VN(N) + IF (N .EQ. 1) RETURN + + DO 1 J=1,N-1 + ED = X-ET(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QN(J) + RETURN + ELSE + QX = QX+QN(J)*C*DY*(X*X-1.D0)/(VN(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INLEGL(N,ET,VN,QN,X,QX) +********************************************************************** +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE LEGENDRE GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + CALL VALEPO(N,X,Y,DY,D2Y) + DN = DFLOAT(N) + C = 1.D0/(DN*(DN+1.D0)) + QX = QN(0)*C*DY*(X-1.D0)/VN(0) + QX = QX+QN(N)*C*DY*(X+1.D0)/VN(N) + IF (N .EQ. 1) RETURN + + DO 1 J=1,N-1 + ED = X-ET(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QN(J) + RETURN + ELSE + QX = QX+QN(J)*C*DY*(X*X-1.D0)/(VN(J)*ED) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE INCHGL(N,QN,X,QX) +********************************************************************* +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE CHEBYSHEV GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + PI = 3.14159265358979323846D0 + CALL VACHPO(N,X,Y,DY,D2Y) + DN = DFLOAT(N) + SN = DFLOAT(1+4*(N/2)-2*N) + PN = PI/DN + C = 1.D0/(DN*DN) + QX = .5D0*SN*QN(0)*C*DY*(X-1.D0) + QX = QX+.5D0*QN(N)*C*DY*(X+1.D0) + IF (N .EQ. 1) RETURN + + SJ = -1.D0 + DO 1 J=1,N-1 + ED = X+DCOS(PN*DFLOAT(J)) + IF(DABS(ED) .LT. EPS) THEN + QX = QN(J) + RETURN + ELSE + QX = QX+QN(J)*SN*SJ*C*DY*(X*X-1.D0)/ED + ENDIF + SJ = -SJ +1 CONTINUE + + RETURN + END +C + SUBROUTINE INLAGR(N,A,ET,VN,QN,X,QX) +********************************************************************* +* COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED +* BY THE VALUES ATTAINED AT THE LAGUERRE GAUSS-RADAU NODES +* N = THE NUMBER OF NODES +* A = PARAMETER > -1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N-1 +* VN = SCALED LAGUERRE FUNCTION AT THE NODES, VN(I), I=0,N-1 +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1 +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* QX = VALUE OF THE POLYNOMIAL IN X +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*) + IF (N .EQ. 0) RETURN + + QX = QN(0) + IF (N .EQ. 1) RETURN + + EPS = 1.D-14 + DN = DFLOAT(N) + C = -1.D0/DN + PR = 1.D0 + SU = 0.D0 + DO 1 M=1,N + DM = 4.D0*DFLOAT(M) + PR = PR*(DM+X)/DM + SU = SU+1.D0/(DM+X) +1 CONTINUE + CALL VALASF(N,A,X,Y,DY) + QX = QN(0)*C*(A+1.D0)*PR*(DY+SU*Y)/VN(0) + DO 2 J=1,N-1 + ED = X-ET(J) + IF(DABS(ED) .LT. EPS) THEN + QX = QN(J) + RETURN + ELSE + PR = 1.D0 + DO 3 K=1,N + DK = 4.D0*DFLOAT(K) + PR = PR*(DK+X)/(DK+ET(J)) +3 CONTINUE + QX = QX+QN(J)*C*PR*(DY+SU*Y)*X/(VN(J)*ED) + ENDIF +2 CONTINUE + + RETURN + END +C + SUBROUTINE NOLEGA(N,QZ,WE,QI,QM) +********************************************************************* +* COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE LEGENDRE ZEROES +* N = THE NUMBER OF ZEROES +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WE = VECTOR OF THE LEGENDRE GAUSS WEIGHTS, WE(I), I=1,N +* QI = INTEGRAL NORM OF THE POLYNOMIAL +* QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE ZEROES +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QZ(1), WE(1) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + SU = 0.D0 + QM = 0.D0 + DO 1 J=1,N + Y = DABS(QZ(J)) + IF(Y .GT. QM) QM=Y + IF(Y .LT. EPS) GOTO 1 + SU = SU+Y*Y*WE(J) +1 CONTINUE + QI = DSQRT(SU) + + RETURN + END +C + SUBROUTINE NOCHGA(N,DZ,QZ,WK,QW,QI,QM) +********************************************************************** +* COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE CHEBYSHEV ZEROES +* N = THE NUMBER OF ZEROES +* DZ = CHEBYSHEV POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WK = VECTOR OF THE CLENSHAW-CURTIS WEIGHTS, WE(I), I=0,2*N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +* QI = INTEGRAL NORM OF THE POLYNOMIAL +* QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE ZEROES +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION DZ(1), QZ(1), WK(0:*) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + DN = DFLOAT(N) + X = -1.D0 + CALL INCHGA(N,DZ,QZ,X,QX) + S1 = 0.D0 + S2 = QX*QX*WK(0) + QM = 0.D0 + DO 1 J=1,N + DJ = DFLOAT(J) + J2 = 2*J + X = -DCOS(PI*DJ/DN) + Y = DABS(QZ(J)) + IF(Y .GT. QM) QM=Y + CALL INCHGA(N,DZ,QZ,X,QX) + S1 = S1+Y*Y + S2 = S2+Y*Y*WK(J2-1)+QX*QX*WK(J2) +1 CONTINUE + QW = DSQRT(S1*PI/DN) + QI = DSQRT(S2) + + RETURN + END +C + SUBROUTINE NOJAGL(N,A,B,VN,QN,WT,QW,QS,QM) +********************************************************************** +* COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE JACOBI GAUSS- +* LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER > -1 +* B = PARAMETER > -1 +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* WT = VECTOR OF THE JACOBI GAUSS-LOBATTO WEIGHTS, WT(I), I=0,N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +* QS = QUADRATURE NORM OF THE POLYNOMIAL +* QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION VN(0:*), QN(0:*), WT(0:*) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + AB2 = AB+2.D0 + DN = DFLOAT(N) + C = ((2.D0)**(AB+1.D0))*(DN+AB+1.D0)/(2.D0*DN+AB+1.D0) + CALL GAMMAF(A1,GA1) + CALL GAMMAF(B1,GB1) + CALL GAMMAF(AB2,GAB2) + C = C*GA1*GB1/GAB2 + DO 1 K=1,N + DK = DFLOAT(K) + C = C*(DK+A)*(DK+B)/(DK*(DK+AB+1.D0)) +1 CONTINUE + + S1 = 0.D0 + S2 = 0.D0 + S3 = 0.D0 + QM = 0.D0 + DO 2 J=0,N + Y1 = QN(J) + YM = DABS(Y1) + IF(YM .GT. QM) QM=YM + Y2 = VN(J) + S2 = S2+Y1*Y2*WT(J) + S3 = S3+Y2*Y2*WT(J) + IF(YM .LT. EPS) GOTO 2 + S1 = S1+Y1*Y1*WT(J) +2 CONTINUE + QS = DSQRT(S1) + QW = DSQRT(S1-(S3-C)*S2*S2/(S3*S3)) + + RETURN + END +C + SUBROUTINE NOLEGL(N,VN,QN,WT,QI,QS,QM) +********************************************************************** +* COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE LEGENDRE GAUSS- +* LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* WT = VECTOR OF THE LEGENDRE GAUSS-LOBATTO WEIGHTS, WT(I), I=0,N +* QW = INTEGRAL NORM OF THE POLYNOMIAL +* QS = QUADRATURE NORM OF THE POLYNOMIAL +* QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION VN(0:*), QN(0:*), WT(0:*) + IF (N .EQ. 0) RETURN + + EPS = 1.D-14 + DN = DFLOAT(N) + C = .5D0*DN*(DN+1.D0)/(2.D0*DN+1.D0) + + S1 = 0.D0 + S2 = 0.D0 + QM = 0.D0 + DO 1 J=0,N + Y1 = QN(J) + YM = DABS(Y1) + IF(YM .GT. QM) QM=YM + Y2 = VN(J) + S2 = S2+Y1*Y2*WT(J) + IF(YM .LT. EPS) GOTO 1 + S1 = S1+Y1*Y1*WT(J) +1 CONTINUE + QS = DSQRT(S1) + QI = DSQRT(DABS(S1-C*S2*S2)) + + RETURN + END +C + SUBROUTINE NOCHGL(N,QN,WK,QW,QI,QS,QM) +********************************************************************** +* COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE CHEBYSHEV GAUSS- +* LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* WK = VECTOR OF THE CLENSHAW-CURTIS WEIGHTS, WE(I), I=0,2*N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +* QI = INTEGRAL NORM OF THE POLYNOMIAL +* QS = QUADRATURE NORM OF THE POLYNOMIAL +* QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*), WK(0:*) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + PH = 1.57079632679489661923D0 + DN = DFLOAT(N) + SN = DFLOAT(1+4*(N/2)-2*N) + Y = QN(0) + S1 = .5D0*Y*Y + S2 = .5D0*Y*SN + S3 = Y*Y*WK(0) + QM = DABS(Y) + + SJ = -1.D0 + DO 1 J=1,N-1 + J2 = 2*J + DJ = DFLOAT(J2-1) + X = -DCOS(PH*DJ/DN) + Y = QN(J) + YM = DABS(Y) + IF(YM .GT. QM) QM=YM + CALL INCHGL(N,QN,X,QX) + S1 = S1+Y*Y + S2 = S2+Y*SN*SJ + S3 = S3+QX*QX*WK(J2-1)+Y*Y*WK(J2) + SJ = -SJ +1 CONTINUE + N2 = 2*N + DD = DFLOAT(N2-1) + X = -DCOS(PH*DD/DN) + Y = QN(N) + YM = DABS(Y) + IF(YM .GT. QM) QM=YM + CALL INCHGL(N,QN,X,QX) + S1 = S1+.5D0*Y*Y + S2 = S2+.5D0*Y + S3 = S3+QX*QX*WK(N2-1)+Y*Y*WK(N2) + + QW = DSQRT(DABS(PI*S1/DN-PH*S2*S2/(DN*DN))) + QI = DSQRT(S3) + QS = DSQRT(PI*S1/DN) + + RETURN + END +C + SUBROUTINE COJAGA(N,A,B,CS,QZ,WE,CO) +********************************************************************** +* COMPUTES THE JACOBI FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE JACOBI ZEROES +* N = THE NUMBER OF ZEROES +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), WE(1), CO(0:*) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + AB2 = AB+2.D0 + CALL GAMMAF(A1,GA1) + CALL GAMMAF(B1,GB1) + CALL GAMMAF(AB2,GAB2) + C = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2 + + SU = 0.D0 + DO 1 J=1,N + SU = SU+QZ(J)*WE(J) + CO(J-1) = 0.D0 +1 CONTINUE + CO(0) = SU/C + IF (N .EQ. 1) RETURN + + DO 2 J=1,N + X = CS(J) + YP = QZ(J)*WE(J) + Y = .5D0*YP*(AB2*X+A-B) + DO 3 K=1,N-1 + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + CC = 2.D0*DK+AB + C1 = 2.D0*DK*(DK+AB)*(CC-2.D0) + C2 = (CC-1.D0)*(CC-2.D0)*CC + C3 = (CC-1.D0)*(A-B)*AB + C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0) + YM = Y + Y = ((C2*X+C3)*Y-C4*YP)/C1 + YP = YM +3 CONTINUE +2 CONTINUE + + DO 4 K=1,N-1 + DK = DFLOAT(K) + C = C*(DK+A)*(DK+B)/DK + CO(K) = CO(K)*(2.D0*DK+AB+1.D0)/C + C = C/(DK+AB+1.D0) +4 CONTINUE + + RETURN + END +C + SUBROUTINE COLEGA(N,CS,QZ,WE,CO) +********************************************************************** +* COMPUTES THE LEGENDRE FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE LEGENDRE ZEROES +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), WE(1), CO(0:*) + IF (N .EQ. 0) RETURN + + SU = 0.D0 + DO 1 J=1,N + SU = SU+QZ(J)*WE(J) + CO(J-1) = 0.D0 +1 CONTINUE + CO(0) = .5D0*SU + IF (N .EQ. 1) RETURN + + DO 2 J=1,N + X = CS(J) + YP = QZ(J)*WE(J) + Y = X*YP + DO 3 K=1,N-1 + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + C1 = 2.D0*DK-1.D0 + C2 = DK-1.D0 + YM = Y + Y = (C1*X*Y-C2*YP)/DK + YP = YM +3 CONTINUE +2 CONTINUE + + DO 4 K=1,N-1 + CO(K) = .5D0*CO(K)*(2.D0*DFLOAT(K)+1.D0) +4 CONTINUE + + RETURN + END +C + SUBROUTINE COCHGA(N,QZ,CO) +********************************************************************** +* COMPUTES THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE CHEBYSHEV ZEROES +* N = THE NUMBER OF ZEROES +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QZ(1), CO(0:*) + IF (N .EQ. 0) RETURN + + PH = 1.57079632679489661923D0 + DN = DFLOAT(N) + SU = 0.D0 + DO 1 J=1,N + SU = SU+QZ(J) +1 CONTINUE + CO(0) = SU/DN + IF (N .EQ. 1) RETURN + + SK = -2.D0 + DO 2 K=1,N-1 + DK = DFLOAT(K) + SU = 0.D0 + DO 3 J=1,N + DJ = 2.D0*DFLOAT(J)-1.D0 + SU = SU+QZ(J)*DCOS(DK*DJ*PH/DN) +3 CONTINUE + CO(K) = SK*SU/DN + SK = -SK +2 CONTINUE + + RETURN + END +C + SUBROUTINE COLAGA(N,A,CS,QZ,WE,CO) +********************************************************************** +* COMPUTES THE LAGUERRE FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE LAGUERRE ZEROES +* N = THE NUMBER OF ZEROES +* A = PARAMETER >-1 +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), WE(1), CO(0:*) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + CALL GAMMAF(A1,C) + SU = 0.D0 + DO 1 J=1,N + SU = SU+QZ(J)*WE(J) + CO(J-1) = 0.D0 +1 CONTINUE + CO(0) = SU/C + IF (N .EQ. 1) RETURN + + DO 2 J=1,N + X = CS(J) + YP = QZ(J)*WE(J) + Y = (A1-X)*YP + DO 3 K=1,N-1 + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + B1 = (2.D0*DK+A-1.D0-X)/DK + B2 = (DK+A-1.D0)/DK + YM = Y + Y = B1*Y-B2*YP + YP = YM +3 CONTINUE +2 CONTINUE + + DO 4 K=1,N-1 + DK = DFLOAT(K) + C = C*(DK+A)/DK + CO(K) = CO(K)/C +4 CONTINUE + + RETURN + END +C + SUBROUTINE COHEGA(N,CS,QZ,WE,CO) +********************************************************************** +* COMPUTES THE HERMITE FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE HERMITE ZEROES +* N = THE NUMBER OF ZEROES +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), WE(1), CO(0:*) + IF (N .EQ. 0) RETURN + + PR = 1.77245385090551588D0 + SU = 0.D0 + DO 1 J=1,N + SU = SU+QZ(J)*WE(J) + CO(J-1) = 0.D0 +1 CONTINUE + CO(0) = SU/PR + IF (N .EQ. 1) RETURN + + DO 2 J=1,N + X = CS(J) + YP = QZ(J)*WE(J)/PR + Y = X*YP + DO 3 K=1,N-1 + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + YM = Y + Y = (X*Y-.5D0*YP)/DK + YP = YM +3 CONTINUE +2 CONTINUE + + RETURN + END +C + SUBROUTINE COJAGL(N,A,B,ET,VN,QN,WT,CO) +********************************************************************** +* COMPUTES THE JACOBI FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY ITS VALUES AT THE JACOBI GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*), WT(0:*), CO(0:*) + CO(0) = QN(0) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + AB2 = AB+2.D0 + CO(1) = (QN(1)-QN(0))/AB2 + CO(0) = .5D0*(QN(0)+QN(1)-(A-B)*CO(1)) + IF (N .EQ. 1) RETURN + + CALL GAMMAF(A1,GA1) + CALL GAMMAF(B1,GB1) + CALL GAMMAF(AB2,GAB2) + C = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2 + SU = 0.D0 + DO 1 J=0,N + SU = SU+QN(J)*WT(J) + CO(J) = 0.D0 +1 CONTINUE + CO(0) = SU/C + + CN = 0.D0 + DO 2 J=0,N + X = ET(J) + YP = QN(J)*WT(J) + Y = .5D0*YP*(AB2*X+A-B) + CN = CN+VN(J)*VN(J)*WT(J) + DO 3 K=1,N + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + CC = 2.D0*DK+AB + C1 = 2.D0*DK*(DK+AB)*(CC-2.D0) + C2 = (CC-1.D0)*(CC-2.D0)*CC + C3 = (CC-1.D0)*(A-B)*AB + C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0) + YM = Y + Y = ((C2*X+C3)*Y-C4*YP)/C1 + YP = YM +3 CONTINUE +2 CONTINUE + + DO 4 K=1,N-1 + DK = DFLOAT(K) + C = C*(DK+A)*(DK+B)/DK + CO(K) = CO(K)*(2.D0*DK+AB+1.D0)/C + C = C/(DK+AB+1.D0) +4 CONTINUE + CO(N) = CO(N)/CN + + RETURN + END +C + SUBROUTINE COLEGL(N,ET,QN,WT,CO) +********************************************************************** +* COMPUTES THE LEGENDRE FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY ITS VALUES AT THE LEGENDRE GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), QN(0:*), WT(0:*), CO(0:*) + CO(0) = QN(0) + IF (N .EQ. 0) RETURN + + CO(0) = .5D0*(QN(0)+QN(1)) + CO(1) = .5D0*(QN(1)-QN(0)) + IF (N .EQ. 1) RETURN + + SU = 0.D0 + DO 1 J=0,N + SU = SU+QN(J)*WT(J) + CO(J) = 0.D0 +1 CONTINUE + CO(0) = .5D0*SU + + DO 2 J=0,N + X = ET(J) + YP = QN(J)*WT(J) + Y = X*YP + DO 3 K=1,N + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + C1 = 2.D0*DK-1.D0 + C2 = DK-1.D0 + YM = Y + Y = (C1*X*Y-C2*YP)/DK + YP = YM +3 CONTINUE +2 CONTINUE + + DN = DFLOAT(N) + CO(N) = .5D0*DN*CO(N) + IF (N .EQ. 1) RETURN + + DO 4 K=1,N-1 + CO(K) = .5D0*CO(K)*(2.D0*DFLOAT(K)+1.D0) +4 CONTINUE + + RETURN + END +C + SUBROUTINE COCHGL(N,QN,CO) +********************************************************************** +* COMPUTES THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY ITS VALUES AT THE CHEBYSHEV GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*), CO(0:*) + CO(0) = QN(0) + IF (N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + DN = DFLOAT(N) + DD = DFLOAT(1+4*(N/2)-2*N) + CO(0) = .5D0*(QN(0)+QN(N)) + CO(N) = .5D0*(QN(0)+DD*QN(N)) + IF (N .EQ. 1) RETURN + + S0 = CO(0) + SN = CO(N) + SJ = -1.D0 + DO 1 J=1,N-1 + S0 = S0+QN(J) + SN = SN+QN(J)*SJ + SJ = -SJ +1 CONTINUE + CO(0) = S0/DN + CO(N) = DD*SN/DN + + SK = -1.D0 + DO 2 K=1,N-1 + DK = DFLOAT(K) + SU = .5D0*(QN(0)+QN(N)*SK) + DO 3 J=1,N-1 + DJ = DFLOAT(J) + SU = SU+QN(J)*DCOS(DK*DJ*PI/DN) +3 CONTINUE + CO(K) = 2.D0*SK*SU/DN + SK = -SK +2 CONTINUE + + RETURN + END +C + SUBROUTINE COLAGR(N,A,ET,QN,WT,CO) +********************************************************************** +* COMPUTES THE LAGUERRE FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY ITS VALUES AT THE LAGUERRE GAUSS-RADAU NODES +* N = THE NUMBER OF NODES +* A = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N-1 +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1 +* WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N-1 +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), QN(0:*), WT(0:*), CO(0:*) + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + CALL GAMMAF(A1,C) + SU = 0.D0 + DO 1 J=0,N-1 + SU = SU+QN(J)*WT(J) + CO(J) = 0.D0 +1 CONTINUE + CO(0) = SU/C + IF (N .EQ. 1) RETURN + + DO 2 J=0,N-1 + X = ET(J) + YP = QN(J)*WT(J) + Y = (A1-X)*YP + DO 3 K=1,N-1 + CO(K) = CO(K)+Y + DK = DFLOAT(K+1) + B1 = (2.D0*DK+A-1.D0-X)/DK + B2 = (DK+A-1.D0)/DK + YM = Y + Y = B1*Y-B2*YP + YP = YM +3 CONTINUE +2 CONTINUE + + DO 4 K=1,N-1 + DK = DFLOAT(K) + C = C*(DK+A)/DK + CO(K) = CO(K)/C +4 CONTINUE + + RETURN + END +C + SUBROUTINE PVJAEX(N,A,B,X,CO,Y,DY,D2Y) +********************************************************************** +* COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND +* SECOND DERIVATIVES BY KNOWING THE JACOBI FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + Y = CO(0) + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + AB = A+B + P = .5D0*((AB+2.D0)*X+A-B) + DP = .5D0*(AB+2.D0) + D2P = 0.D0 + Y = CO(0)+P*CO(1) + DY = DP*CO(1) + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + PP = 1.D0 + DPP = 0.D0 + D2PP = 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + CC = 2.D0*DK+AB + C1 = 2.D0*DK*(DK+AB)*(CC-2.D0) + C2 = (CC-1.D0)*(CC-2.D0)*CC + C3 = (CC-1.D0)*(A-B)*AB + C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0) + PM = P + P = ((C2*X+C3)*P-C4*PP)/C1 + Y = Y+P*CO(K) + PP = PM + DPM = DP + DP = ((C2*X+C3)*DP-C4*DPP+C2*PP)/C1 + DY = DY+DP*CO(K) + DPP = DPM + D2PM = D2P + D2P = ((C2*X+C3)*D2P-C4*D2PP+2.D0*C2*DPP)/C1 + D2Y = D2Y+D2P*CO(K) + D2PP = D2PM +1 CONTINUE + + RETURN + END +C + SUBROUTINE PVLEEX(N,X,CO,Y,DY,D2Y) +********************************************************************** +* COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND +* SECOND DERIVATIVES BY KNOWING THE LEGENDRE FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + Y = CO(0) + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + P = X + DP = 1.D0 + D2P = 0.D0 + Y = CO(0)+P*CO(1) + DY = DP*CO(1) + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + PP = 1.D0 + DPP = 0.D0 + D2PP = 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + C2 = 2.D0*DK-1.D0 + C4 = DK-1.D0 + PM = P + P = (C2*X*P-C4*PP)/DK + Y = Y+P*CO(K) + PP = PM + DPM = DP + DP = (C2*X*DP-C4*DPP+C2*PP)/DK + DY = DY+DP*CO(K) + DPP = DPM + D2PM = D2P + D2P = (C2*X*D2P-C4*D2PP+2.D0*C2*DPP)/DK + D2Y = D2Y+D2P*CO(K) + D2PP = D2PM +1 CONTINUE + + RETURN + END +C + SUBROUTINE PVCHEX(N,X,CO,Y,DY,D2Y) +********************************************************************** +* COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND +* SECOND DERIVATIVES BY KNOWING THE CHEBYSHEV FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + Y = CO(0) + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + P = X + DP = 1.D0 + D2P = 0.D0 + Y = CO(0)+P*CO(1) + DY = DP*CO(1) + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + PP = 1.D0 + DPP = 0.D0 + D2PP = 0.D0 + DO 1 K=2,N + PM = P + P = 2.D0*X*P-PP + Y = Y+P*CO(K) + PP = PM + DPM = DP + DP = 2.D0*X*DP+2.D0*PP-DPP + DY = DY+DP*CO(K) + DPP = DPM + D2PM = D2P + D2P = 2.D0*X*D2P+4.D0*DPP-D2PP + D2Y = D2Y+D2P*CO(K) + D2PP = D2PM +1 CONTINUE + + RETURN + END +C + SUBROUTINE PVLAEX(N,A,X,CO,Y,DY,D2Y) +********************************************************************** +* COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND +* SECOND DERIVATIVES BY KNOWING THE LAGUERRE FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + Y = CO(0) + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + P = 1.D0+A-X + DP = -1.D0 + D2P = 0.D0 + Y = CO(0)+P*CO(1) + DY = DP*CO(1) + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + PP = 1.D0 + DPP = 0.D0 + D2PP = 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + B1 = (2.D0*DK+A-1.D0-X)/DK + B2 = (DK+A-1.D0)/DK + PM = P + P = B1*P-B2*PP + Y = Y+P*CO(K) + PP = PM + DPM = DP + DP = B1*DP-PP/DK-B2*DPP + DY = DY+DP*CO(K) + DPP = DPM + D2PM = D2P + D2P = B1*D2P-2.D0*DPP/DK-B2*D2PP + D2Y = D2Y+D2P*CO(K) + D2PP = D2PM +1 CONTINUE + + RETURN + END +C + SUBROUTINE PVHEEX(N,X,CO,Y,DY,D2Y) +********************************************************************** +* COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND +* SECOND DERIVATIVES BY KNOWING THE HERMITE FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* X = THE POINT IN WHICH THE COMPUTATION IS PERFORMED +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* Y = VALUE OF THE POLYNOMIAL IN X +* DY = VALUE OF THE FIRST DERIVATIVE IN X +* D2Y= VALUE OF THE SECOND DERIVATIVE IN X +********************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + Y = CO(0) + DY = 0.D0 + D2Y = 0.D0 + IF (N .EQ. 0) RETURN + + P = 2.D0*X + DP = 2.D0 + D2P = 0.D0 + Y = CO(0)+P*CO(1) + DY = DP*CO(1) + D2Y = 0.D0 + IF (N .EQ. 1) RETURN + + PP = 1.D0 + DPP = 0.D0 + D2PP = 0.D0 + DO 1 K=2,N + DK = DFLOAT(K) + PM = P + P = 2.D0*X*P-2.D0*PP*(DK-1.D0) + Y = Y+P*CO(K) + DY = DY+2.D0*DK*PM*CO(K) + D2Y = D2Y+4.D0*DK*(DK-1.D0)*PP*CO(K) + PP = PM +1 CONTINUE + + RETURN + END +C + SUBROUTINE NOJAEX(N,A,B,CO,QW) +*************************************************************** +* COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE +* FOURIER COEFFICIENTS WITH RESPECT TO THE JACOBI BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +*************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + + EPS = 1.D-14 + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + AB2 = AB+2.D0 + DN = DFLOAT(N) + CALL GAMMAF(A1,GA1) + CALL GAMMAF(B1,GB1) + CALL GAMMAF(AB2,GAB2) + C = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2 + V = DABS(CO(0)) + QW = V*DSQRT(C) + IF (N .EQ. 0) RETURN + + SU = 0.D0 + IF (V .LT. EPS) GOTO 1 + SU = C*V*V +1 DO 2 K=1,N + DK = DFLOAT(K) + C = C*(DK+A)*(DK+B)/DK + V = DABS(CO(K)) + IF (V .LT. EPS) GOTO 3 + SU = SU+C*V*V/(2.D0*DK+AB+1.D0) +3 C = C/(DK+AB+1.D0) +2 CONTINUE + QW = DSQRT(SU) + + RETURN + END +C + SUBROUTINE NOLEEX(N,CO,QI) +**************************************************************** +* COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE +* FOURIER COEFFICIENTS WITH RESPECT TO THE LEGENDRE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QI = INTEGRAL NORM OF THE POLYNOMIAL +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + + EPS = 1.D-14 + SU = 0.D0 + DO 1 K=0,N + DK = DFLOAT(K) + V = DABS(CO(K)) + IF (V .LT. EPS) GOTO 1 + SU = SU+V*V/(2.D0*DK+1.D0) +1 CONTINUE + QI = DSQRT(2.D0*SU) + + RETURN + END +C + SUBROUTINE NOCHEX(N,CO,QW,QI) +**************************************************************** +* COMPUTES THE INTEGRAL NORMS OF A POLYNOMIAL BY KNOWING THE +* FOURIER COEFFICIENTS WITH RESPECT TO THE CHEBYSHEV BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +* QI = INTEGRAL NORM OF THE POLYNOMIAL +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + + EPS = 1.D-14 + PR = 1.77245385090551588D0 + R2 = 1.41421356237309515D0 + V = DABS(CO(0)) + QW = PR*V + QI = R2*V + IF (N .EQ. 0) RETURN + + SU = 0.D0 + IF (V .LT. EPS) GOTO 1 + SU = 2.D0*V*V +1 DO 2 K=1,N + V = DABS(CO(K)) + IF (V .LT. EPS) GOTO 2 + SU = SU+V*V +2 CONTINUE + QW = PR*DSQRT(.5D0*SU) + + SU = 0.D0 + DO 3 K=0,N,2 + V = CO(K) + DO 3 M=0,N,2 + D1 = 1.D0-DFLOAT((K+M)*(K+M)) + D2 = 1.D0-DFLOAT((K-M)*(K-M)) + C =1.D0/D1+1.D0/D2 + SU = SU+C*V*CO(M) +3 CONTINUE + DO 4 K=1,N,2 + V = CO(K) + DO 4 M=1,N,2 + D1 = 1.D0-DFLOAT((K+M)*(K+M)) + D2 = 1.D0-DFLOAT((K-M)*(K-M)) + C =1.D0/D1+1.D0/D2 + SU = SU+C*V*CO(M) +4 CONTINUE + QI = DSQRT(SU) + + RETURN + END +C + SUBROUTINE NOLAEX(N,A,CO,QW) +***************************************************************** +* COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE +* FOURIER COEFFICIENTS WITH RESPECT TO THE LAGUERRE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +***************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + + EPS = 1.D-14 + A1 = A+1.D0 + CALL GAMMAF(A1,C) + V = DABS(CO(0)) + QW = V*DSQRT(C) + IF (N .EQ. 0) RETURN + + SU = 0.D0 + IF (V .LT. EPS) GOTO 1 + SU = C*V*V +1 DO 2 K=1,N + DK = DFLOAT(K) + C = C*(DK+A)/DK + V = DABS(CO(K)) + IF (V .LT. EPS) GOTO 2 + SU = SU+C*V*V +2 CONTINUE + QW = DSQRT(SU) + + RETURN + END +C + SUBROUTINE NOHEEX(N,CO,QW) +**************************************************************** +* COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE +* FOURIER COEFFICIENTS WITH RESPECT TO THE HERMITE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL +**************************************************************** + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*) + + EPS = 1.D-14 + PR = 1.33133536380038953D0 + R2 = 1.41421356237309515D0 + V = DABS(CO(0)) + QW = V*PR + IF (N .EQ. 0) RETURN + + SU = 0.D0 + IF (V .LT. EPS) GOTO 1 + SU = V*V +1 C = 1.D0 + DO 2 K=1,N + DK = DFLOAT(K) + C = C*R2*DSQRT(DK) + V = DABS(CO(K)) + IF (V .LT. EPS) GOTO 2 + SU = SU+(C*V)*(C*V) +2 CONTINUE + QW = PR*DSQRT(SU) + + RETURN + END +C + SUBROUTINE COJADE(N,G,CO,CD,CD2) +************************************************************************ +* COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL +* FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE JACOBI BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* G = PARAMETER >-1 +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* CD = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N +* CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*), CD(0:*), CD2(0:*) + + CD(N) = 0.D0 + CD2(N) = 0.D0 + IF (N .EQ. 0) RETURN + + CD(0) = (G+1.D0)*CO(1) + CD2(N-1) = 0.D0 + IF (N .EQ. 1) RETURN + + DN = DFLOAT(N) + G2 = 2.D0*G + CD(N-1) = (2.D0*DN+G2-1.D0)*(DN+G)*CO(N)/(DN+G2) + DO 1 K=0,N-2 + KR = N-K-2 + IF (KR .NE. 0) THEN + DK = DFLOAT(KR) + C1 = (2.D0*DK+G2+1.D0)*(DK+G+1.D0)/(DK+G2+1.D0) + C2 = (DK+G+2.D0)/((2.D0*DK+G2+5.D0)*(DK+G2+2.D0)) + CD(KR) = C1*(C2*CD(KR+2)+CO(KR+1)) + CD2(KR) = C1*(C2*CD2(KR+2)+CD(KR+1)) + ELSE + CD(0) = .25D0*(G+2.D0)*CD(2)/(G+2.5D0)+(G+1.D0)*CO(1) + CD2(0) = .25D0*(G+2.D0)*CD2(2)/(G+2.5D0)+(G+1.D0)*CD(1) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE COLEDE(N,CO,CD,CD2) +************************************************************************ +* COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL +* FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE LEGENDRE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* CD = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N +* CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*), CD(0:*), CD2(0:*) + + CD(N) = 0.D0 + CD2(N) = 0.D0 + IF (N .EQ. 0) RETURN + + DN = DFLOAT(N) + CD(N-1) = (2.D0*DN-1.D0)*CO(N) + CD2(N-1) = 0.D0 + IF (N .EQ. 1) RETURN + + DO 1 K=0,N-2 + KR = N-K-2 + DK = 2.D0*DFLOAT(KR)+1.D0 + CD(KR) = DK*(CD(KR+2)/(DK+4.D0)+CO(KR+1)) + CD2(KR) = DK*(CD2(KR+2)/(DK+4.D0)+CD(KR+1)) +1 CONTINUE + + RETURN + END +C + SUBROUTINE COCHDE(N,CO,CD,CD2) +************************************************************************ +* COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL +* FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE CHEBYSHEV BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* CD = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N +* CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*), CD(0:*), CD2(0:*) + + CD(N) = 0.D0 + CD2(N) = 0.D0 + IF (N .EQ. 0) RETURN + + CD(0) = CO(1) + CD2(N-1) = 0.D0 + IF (N .EQ. 1) RETURN + + DN = DFLOAT(N) + CD(N-1) = 2.D0*DN*CO(N) + DO 1 K=0,N-2 + KR = N-K-2 + IF (KR .NE. 0) THEN + DK = 2.D0*(DFLOAT(KR)+1.D0) + CD(KR) = CD(KR+2)+DK*CO(KR+1) + CD2(KR) = CD2(KR+2)+DK*CD(KR+1) + ELSE + CD(0) = .5D0*CD(2)+CO(1) + CD2(0) = .5D0*CD2(2)+CD(1) + ENDIF +1 CONTINUE + + RETURN + END +C + SUBROUTINE COLADE(N,CO,CD,CD2) +************************************************************************ +* COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL +* FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE LAGUERRE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* CD = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N +* CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*), CD(0:*), CD2(0:*) + + CD(N) = 0.D0 + CD2(N) = 0.D0 + IF (N .EQ. 0) RETURN + + CD(N-1) = -CO(N) + CD2(N-1) = 0.D0 + IF (N .EQ. 1) RETURN + + DO 1 K=0,N-2 + KR = N-K-2 + CD(KR) = CD(KR+1)-CO(KR+1) + CD2(KR) = CD2(KR+2)-CD(KR+1) +1 CONTINUE + + RETURN + END +C + SUBROUTINE COHEDE(N,CO,CD,CD2) +************************************************************************ +* COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL +* FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE HERMITE BASIS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* CD = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N +* CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CO(0:*), CD(0:*), CD2(0:*) + + CD(N) = 0.D0 + CD2(N) = 0.D0 + IF (N .EQ. 0) RETURN + + DN = DFLOAT(N) + CD(N-1) = 2.D0*DN*CO(N) + CD2(N-1) = 0.D0 + IF (N .EQ. 1) RETURN + + DO 1 K=0,N-2 + KR = N-K-2 + DK = 2.D0*DFLOAT(KR)+2.D0 + CD(KR) = DK*CO(KR+1) + CD2(KR) = DK*CD(KR+1) +1 CONTINUE + + RETURN + END +C + SUBROUTINE DEJAGA(N,A,B,CS,DZ,QZ,DQZ) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE JACOBI ZEROES FROM +* THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE NUMBER OF ZEROES +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* CS = ZEROES OF THE JACOBI POLYNOMIAL, CS(I), I=1,N +* DZ = JACOBI DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), QZ(1), DQZ(1) + IF (N .EQ. 0) RETURN + + DO 1 I=1,N + SU = 0.D0 + CI = CS(I) + DI = DZ(I) + DO 2 J=1,N + IF (I .NE. J) THEN + CJ = CS(J) + DJ = DZ(J) + SU = SU+QZ(J)/(DJ*(CI-CJ)) + ELSE + SU = SU+.5D0*QZ(I)*((A+B+2.D0)*CI+A-B)/(DI*(1.D0-CI*CI)) + ENDIF +2 CONTINUE + DQZ(I) = DI*SU +1 CONTINUE + + RETURN + END +C + SUBROUTINE DELAGA(N,A,CS,QZ,DQZ) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LAGUERRE ZEROES FROM +* THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE NUMBER OF ZEROES +* A = PARAMETER >-1 +* CS = ZEROES OF THE LAGUERRE POLYNOMIAL, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), DQZ(1) + IF (N .EQ. 0) RETURN + + DO 1 J=1,N + CJ = CS(J) + CALL VALAPO(N,A,CJ,Y,DY,D2Y) + QZ(J) = QZ(J)/DY +1 CONTINUE + + DO 2 I=1,N + SU = 0.D0 + CI = CS(I) + CALL VALAPO(N,A,CI,Y,DI,D2Y) + DO 3 J=1,N + IF (I .NE. J) THEN + CJ = CS(J) + SU = SU+QZ(J)/(CI-CJ) + ELSE + SU = SU+.5D0*QZ(I)*(CI-A-1.D0)/CI + ENDIF +3 CONTINUE + DQZ(I) = DI*SU +2 CONTINUE + + DO 4 I=1,N + CI = CS(I) + CALL VALAPO(N,A,CI,Y,DY,D2Y) + QZ(I) = DY*QZ(I) +4 CONTINUE + + RETURN + END +C + SUBROUTINE DEHEGA(N,CS,QZ,DQZ) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE HERMITE ZEROES FROM +* THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE NUMBER OF ZEROES +* CS = ZEROES OF THE HERMITE POLYNOMIAL, CS(I), I=1,N +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), QZ(1), DQZ(1) + IF (N .EQ. 0) RETURN + + DO 1 J=1,N + CJ = CS(J) + CALL VAHEPO(N,CJ,Y,DY,D2Y) + QZ(J) = QZ(J)/DY +1 CONTINUE + + DO 2 I=1,N + SU = 0.D0 + CI = CS(I) + CALL VAHEPO(N,CI,Y,DI,D2Y) + DO 3 J=1,N + IF (I .NE. J) THEN + CJ = CS(J) + SU = SU+QZ(J)/(CI-CJ) + ELSE + SU = SU+CI*QZ(I) + ENDIF +3 CONTINUE + DQZ(I) = DI*SU +2 CONTINUE + + DO 4 I=1,N + CI = CS(I) + CALL VAHEPO(N,CI,Y,DY,D2Y) + QZ(I) = DY*QZ(I) +4 CONTINUE + + RETURN + END +C + SUBROUTINE DEJAGL(N,A,B,ET,VN,QN,DQN) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE JACOBI GAUSS-LOBATTO +* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE DEGREE OF THE POLYNOMIAL +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*) + DQN(0) = 0.D0 + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + DN = DFLOAT(N) + C1 = DN*(DN+AB+1.D0) + C2 = A1*VN(0)/(B1*VN(N)) + DQN(0) = .5D0*((A-C1)*QN(0)/(B+2.D0)-C2*QN(N)) + DQN(N) = .5D0*(QN(0)/C2+(C1-B)*QN(N)/(A+2.D0)) + IF (N .EQ. 1) RETURN + + S1 = DQN(0) + S2 = DQN(N) + C3 = -VN(0)/B1 + C4 = VN(N)/A1 + DO 1 J=1,N-1 + VJ = QN(J)/VN(J) + EJ = ET(J) + S1 = S1+C3*VJ/(1.D0+EJ) + S2 = S2+C4*VJ/(1.D0-EJ) +1 CONTINUE + DQN(0) = S1 + DQN(N) = S2 + + DO 2 I=1,N-1 + VI = VN(I) + EI = ET(I) + SU = B1*QN(0)/((1.D0+EI)*VN(0))-A1*QN(N)/((1.D0-EI)*VN(N)) + DO 3 J=1,N-1 + IF (I .NE. J) THEN + VJ = VN(J) + EJ = ET(J) + SU = SU+QN(J)/(VJ*(EI-EJ)) + ELSE + SU = SU+.5D0*QN(I)*(AB*EI+A-B)/(VI*(1.D0-EI*EI)) + ENDIF +3 CONTINUE + DQN(I) = VI*SU +2 CONTINUE + + RETURN + END +C + SUBROUTINE DELEGL(N,ET,VN,QN,DQN) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LEGENDRE GAUSS-LOBATTO +* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE DEGREE OF THE POLYNOMIAL +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*) + DQN(0) = 0.D0 + IF (N .EQ. 0) RETURN + + DO 1 I=0,N + SU = 0.D0 + VI = VN(I) + EI = ET(I) + DO 2 J=0,N + IF (I .EQ. J) GOTO 2 + VJ = VN(J) + EJ = ET(J) + SU = SU+QN(J)/(VJ*(EI-EJ)) +2 CONTINUE + DQN(I) = VI*SU +1 CONTINUE + + DN = DFLOAT(N) + C = .25D0*DN*(DN+1.D0) + DQN(0) = DQN(0)-C*QN(0) + DQN(N) = DQN(N)+C*QN(N) + + RETURN + END +C + SUBROUTINE DECHGL(N,ET,QN,DQN) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE CHEBYSHEV GAUSS-LOBATTO +* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE DEGREE OF THE POLYNOMIAL +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), QN(0:*), DQN(0:*) + DQN(0) = 0.D0 + IF (N .EQ. 0) RETURN + + DN = DFLOAT(N) + CN = (2.D0*DN*DN+1.D0)/6.D0 + SN = DFLOAT(1+4*(N/2)-2*N) + DQN(0) = -CN*QN(0)-.5D0*SN*QN(N) + DQN(N) = .5D0*SN*QN(0)+CN*QN(N) + IF (N .EQ. 1) RETURN + + S1 = DQN(0) + S2 = DQN(N) + SGN = -1.D0 + DO 1 J=1,N-1 + EJ = ET(J) + QJ = 2.D0*SGN*QN(J) + S1 = S1-QJ/(1.D0+EJ) + S2 = S2+QJ*SN/(1.D0-EJ) + SGN = -SGN +1 CONTINUE + DQN(0) = S1 + DQN(N) = S2 + + SGNI = -1.D0 + DO 2 I=1,N-1 + EI = ET(I) + SU = .5D0*SGNI*(QN(0)/(1.D0+EI)-SN*QN(N)/(1.D0-EI)) + SGNJ = -1.D0 + DO 3 J=1,N-1 + IF (I .NE. J) THEN + EJ = ET(J) + SU = SU+SGNI*SGNJ*QN(J)/(EI-EJ) + ELSE + SU = SU-.5D0*EI*QN(I)/(1.D0-EI*EI) + ENDIF + SGNJ = -SGNJ +3 CONTINUE + DQN(I) = SU + SGNI = -SGNI +2 CONTINUE + + RETURN + END +C + SUBROUTINE DELAGR(N,A,ET,QN,DQN) +************************************************************************ +* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LAGUERRE GAUSS-RADAU +* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS +* N = THE NUMBER OF NODES +* A = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N-1 +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1 +* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N-1 +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), QN(0:*), DQN(0:*) + DN = DFLOAT(N) + DQN(0) = (1.D0-DN)*QN(0)/(A+2.D0) + IF (N .EQ. 1) RETURN + + A1 = A+1.D0 + SU = DQN(0) + X = 0.D0 + CALL VALAPO(N,A,X,Y,DY,D2Y) + C = Y + DO 1 J=1,N-1 + EJ = ET(J) + CALL VALAPO(N,A,EJ,Y,DY,D2Y) + QN(J) = QN(J)/Y + SU = SU-C*QN(J)/(A1*EJ) +1 CONTINUE + DQN(0) = SU + + DO 2 I=1,N-1 + EI = ET(I) + CALL VALAPO(N,A,EI,Y,DY,D2Y) + SU = A1*QN(0)/(C*EI) + DO 3 J=1,N-1 + IF (I .NE. J) THEN + EJ = ET(J) + SU = SU+QN(J)/(EI-EJ) + ELSE + SU = SU+.5D0*QN(I)*(EI-A)/EI + ENDIF +3 CONTINUE + DQN(I) = Y*SU +2 CONTINUE + + DO 4 J=1,N-1 + EJ = ET(J) + CALL VALAPO(N,A,EJ,Y,DY,D2Y) + QN(J) = Y*QN(J) +4 CONTINUE + + RETURN + END +C + SUBROUTINE DMJAGL(N,NM,A,B,ET,VN,DMA) +************************************************************************ +* COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE +* JACOBI GAUSS-LOBATTO NODES +* N = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX +* NM = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT +* A = PARAMETER >-1 +* B = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N +* DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N J=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), DMA(0:NM,0:*) + DMA(0,0) = 0.D0 + IF (N .EQ. 0) RETURN + + A1 = A+1.D0 + B1 = B+1.D0 + AB = A+B + DN = DFLOAT(N) + C1 = DN*(DN+AB+1.D0) + C2 = A1*VN(0)/(B1*VN(N)) + DMA(0,0) = .5D0*(A-C1)/(B+2.D0) + DMA(N,N) = .5D0*(C1-B)/(A+2.D0) + DMA(0,N) = -.5D0*C2 + DMA(N,0) = .5D0/C2 + IF (N .EQ. 1) RETURN + + C3 = VN(0)/B1 + C4 = VN(N)/A1 + DO 1 J=1,N-1 + VJ = VN(J) + EJ = ET(J) + DMA(0,J) = -C3/(VJ*(1.D0+EJ)) + DMA(N,J) = C4/(VJ*(1.D0-EJ)) + DMA(J,0) = VJ/(C3*(1.D0+EJ)) + DMA(J,N) = -VJ/(C4*(1.D0-EJ)) +1 CONTINUE + + DO 2 I=1,N-1 + VI = VN(I) + EI = ET(I) + DO 3 J=1,N-1 + IF (I .NE. J) THEN + VJ = VN(J) + EJ = ET(J) + DMA(I,J) = VI/(VJ*(EI-EJ)) + ELSE + DMA(I,I) = .5D0*(AB*EI+A-B)/(1.D0-EI*EI) + ENDIF +3 CONTINUE +2 CONTINUE + + RETURN + END +C + SUBROUTINE DMLEGL(N,NM,ET,VN,DMA) +************************************************************************ +* COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE +* LEGENDRE GAUSS-LOBATTO NODES +* N = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX +* NM = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N J=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), DMA(0:NM,0:*) + DMA(0,0) = 0.D0 + IF (N .EQ. 0) RETURN + + DO 1 I=0,N + VI = VN(I) + EI = ET(I) + DO 2 J=0,N + IF (I .NE. J) THEN + VJ = VN(J) + EJ = ET(J) + DMA(I,J) = VI/(VJ*(EI-EJ)) + ELSE + DMA(I,I) = 0.D0 + ENDIF +2 CONTINUE +1 CONTINUE + + DN = DFLOAT(N) + C = .25D0*DN*(DN+1.D0) + DMA(0,0) = -C + DMA(N,N) = C + + RETURN + END +C + SUBROUTINE DMCHGL(N,NM,ET,DMA) +************************************************************************ +* COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE +* CHEBYSHEV GAUSS-LOBATTO NODES +* N = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX +* NM = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N J=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), DMA(0:NM,0:*) + DMA(0,0) = 0.D0 + IF (N .EQ. 0) RETURN + + DN = DFLOAT(N) + CN = (2.D0*DN*DN+1.D0)/6.D0 + SN = DFLOAT(1+4*(N/2)-2*N) + DMA(0,0) = -CN + DMA(N,N) = CN + DMA(0,N) = -.5D0*SN + DMA(N,0) = .5D0*SN + IF (N .EQ. 1) RETURN + + SGN = -1.D0 + DO 1 J=1,N-1 + EJ = ET(J) + DMA(0,J) = -2.D0*SGN/(1.D0+EJ) + DMA(N,J) = 2.D0*SGN*SN/(1.D0-EJ) + DMA(J,0) = .5D0*SGN/(1.D0+EJ) + DMA(J,N) = -.5D0*SGN*SN/(1.D0-EJ) + SGN = -SGN +1 CONTINUE + + SGNI = -1.D0 + DO 2 I=1,N-1 + EI = ET(I) + SGNJ = -1.D0 + DO 3 J=1,N-1 + IF (I .NE. J) THEN + EJ = ET(J) + DMA(I,J) = SGNI*SGNJ/(EI-EJ) + ELSE + DMA(I,I) = -.5D0*EI/(1.D0-EI*EI) + ENDIF + SGNJ = -SGNJ +3 CONTINUE + SGNI = -SGNI +2 CONTINUE + + RETURN + END +C + SUBROUTINE DMLAGR(N,NM,A,ET,DMA) +************************************************************************ +* COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE +* LAGUERRE GAUSS-RADAU NODES +* N = DIMENSION OF THE MATRIX +* NM = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT +* A = PARAMETER >-1 +* ET = VECTOR OF THE NODES, ET(I), I=0,N-1 +* DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N-1 J=0,N-1 +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), DMA(0:NM,0:*) + DN = DFLOAT(N) + DMA(0,0) = (1.D0-DN)/(A+2.D0) + IF (N .EQ. 1) RETURN + + A1 = A+1.D0 + X = 0.D0 + CALL VALAPO(N,A,X,Y,DY,D2Y) + C = Y + DO 1 J=1,N-1 + EJ = ET(J) + CALL VALAPO(N,A,EJ,Y,DY,D2Y) + DMA(0,J) = -C/(A1*EJ*Y) + DMA(J,0) = A1*Y/(C*EJ) +1 CONTINUE + + DO 2 I=1,N-1 + EI = ET(I) + CALL VALAPO(N,A,EI,Y,DY,D2Y) + DO 3 J=1,N-1 + IF (I .NE. J) THEN + EJ = ET(J) + DMA(I,J) = Y/(EI-EJ) + ELSE + DMA(I,I) = .5D0*(EI-A)/EI + ENDIF +3 CONTINUE +2 CONTINUE + + DO 4 J=1,N-1 + EJ = ET(J) + CALL VALAPO(N,A,EJ,Y,DY,D2Y) + DO 5 I=1,N-1 + IF (I .EQ. J) GOTO 5 + DMA(I,J) = DMA(I,J)/Y +5 CONTINUE +4 CONTINUE + + RETURN + END +C + SUBROUTINE FCCHGA(N,QZ,CO) +************************************************************************ +* COMPUTES USING FFT THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY THE VALUES ATTAINED AT THE CHEBYSHEV ZEROES +* N = THE NUMBER OF ZEROES +* QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1 +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QZ(1), CO(0:*) + IF(N .EQ. 0) RETURN + + PH = 1.57079632679489661923D0 + R2 = 1.41421356237309515D0 + DN = DFLOAT(N) + CO(0) = QZ(1)/DN + IF(N .EQ. 1) RETURN + + N2 = N/2 + C = 2.D0/DSQRT(DN) + SN = DFLOAT(1+4*N2-2*N) + CO(N-N2-1) = QZ(N) + DO 1 I=1,N2 + CO(I-1) = QZ(2*I-1) + CO(N-I) = QZ(2*I) +1 CONTINUE + + CALL C06EAF(CO,N,IFAIL) + IF (IFAIL .NE. 0) THEN + WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EAF' + ENDIF + + CO(0) = .5D0*C*CO(0) + IF (2*N2 .EQ. N) THEN + CO(N2) = C*((-1.D0)**N2)*CO(N2)/R2 + ENDIF + IF (N .EQ. 2) RETURN + SM = -1.D0 + DO 2 M=1,N-N2-1 + AR = PH*DFLOAT(M)/DN + CS = DCOS(AR) + SI = DSIN(AR) + V1 = C*SM*(CO(M)*CS+CO(N-M)*SI) + V2 = C*SM*SN*(CO(M)*SI-CO(N-M)*CS) + CO(M) = V1 + CO(N-M) = V2 + SM = -SM +2 CONTINUE + + RETURN + END +C + SUBROUTINE FCCHGL(N,QN,CO) +************************************************************************ +* COMPUTES USING FFT THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL +* INDIVIDUATED BY ITS VALUES AT THE CHEBYSHEV GAUSS-LOBATTO NODES +* N = THE DEGREE OF THE POLYNOMIAL +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*), CO(0:*) + IF(N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + DN = DFLOAT(N) + N2 = N/2 + SN = DFLOAT(1+4*N2-2*N) + S1 = .5D0*(QN(0)+QN(N)) + S2 = .5D0*(SN*QN(0)+QN(N)) + CO(0) = S1 + CO(N) = S2 + IF(N .EQ. 1) RETURN + + CO(0) = QN(0) + IF (2*N2 .EQ. N) THEN + CO(N2) = QN(N) + ENDIF + IF(N .EQ. 2) GOTO 2 + + DO 1 I=1,N-N2-1 + I2 = 2*I + CO(I) = QN(I2) + CO(N-I) = QN(I2-1)-QN(I2+1) +1 CONTINUE + +2 CALL C06EBF(CO,N,IFAIL) + IF (IFAIL .NE. 0) THEN + WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EBF' + ENDIF + + SJ = -1.D0 + DO 3 J=1,N-1 + S1 = S1+QN(J) + S2 = S2+SJ*SN*QN(J) + SJ = -SJ +3 CONTINUE + CO(0) = S1/DN + CO(N) = S2/DN + C = .5D0/DSQRT(DN) + IF (2*N2 .EQ. N) THEN + CO(N2) = 2.D0*C*((-1.D0)**N2)*CO(N2) + ENDIF + IF(N .EQ. 2) RETURN + + SM = -1.D0 + DO 4 M=1,N-N2-1 + AR = PI*DFLOAT(M)/DN + SI = .5D0/DSIN(AR) + V1 = CO(M)*(1.D0+SI)+CO(N-M)*(1.D0-SI) + V2 = CO(M)*(1.D0-SI)+CO(N-M)*(1.D0+SI) + CO(M) = C*SM*V1 + CO(N-M) = C*SM*SN*V2 + SM = -SM +4 CONTINUE + + RETURN + END +C + SUBROUTINE FVCHGL(N,CO,QN) +******************************************************************* +* COMPUTES USING FFT THE VALUES OF A POLYNOMIAL AT THE CHEBYSHEV +* GAUSS-LOBATTO NODES FROM ITS CHEBYSHEV FOURIER COEFFICIENTS +* N = THE DEGREE OF THE POLYNOMIAL +* CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +******************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*), CO(0:*) + IF(N .EQ. 0) RETURN + + PI = 3.14159265358979323846D0 + DN = DFLOAT(N) + N2 = N/2 + SN = DFLOAT(1+4*N2-2*N) + S1 = CO(0)+SN*CO(N) + S2 = CO(0)+CO(N) + QN(0) = S1 + QN(N) = S2 + IF(N .EQ. 1) RETURN + + QN(0) = 2.D0*CO(0) + IF (2*N2 .EQ. N) THEN + QN(N2) = 2.D0*CO(N) + ENDIF + IF(N .EQ. 2) GOTO 2 + + DO 1 I=1,N-N2-1 + I2 = 2*I + QN(N-I) = CO(I2+1)-CO(I2-1) + QN(I) = CO(I2) +1 CONTINUE + IF (2*N2 .NE. N) THEN + QN(N2+1) = 2.D0*CO(N)-CO(N-2) + ENDIF + +2 CALL C06EBF(QN,N,IFAIL) + IF (IFAIL .NE. 0) THEN + WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EBF' + ENDIF + + SJ = -1.D0 + DO 3 J=1,N-1 + S1 = S1+SJ*CO(J) + S2 = S2+CO(J) + SJ = -SJ +3 CONTINUE + QN(0) = S1 + QN(N) = S2 + C = .25D0*DSQRT(DN) + IF (2*N2 .EQ. N) THEN + QN(N2) = 2.D0*C*QN(N2) + ENDIF + IF(N .EQ. 2) RETURN + + DO 4 M=1,N-N2-1 + AR = PI*DFLOAT(M)/DN + SI = .5D0/DSIN(AR) + V1 = QN(M)*(1.D0+SI)+QN(N-M)*(1.D0-SI) + V2 = QN(M)*(1.D0-SI)+QN(N-M)*(1.D0+SI) + QN(M) = C*V1 + QN(N-M) = C*V2 +4 CONTINUE + + RETURN + END +C + SUBROUTINE FDCHGL(N,QN,CD,DQN) +********************************************************************* +* COMPUTES USING FFT THE FOURIER COEFFICIENTS AND THE VALUES AT THE +* CHEBYSHEV GAUSS-LOBATTO NODES OF THE DERIVATIVE OF A POLYNOMIAL +* FROM THE VALUES ATTAINED BY THE POLYNOMIAL AT THE NODES +* N = THE DEGREE OF THE POLYNOMIAL +* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N +* CD = FOURIER COEFFICIENTS OF THE DERIVATIVE, CD(I), I=0,N +* DQN = VALUES OF THE DERIVATIVE AT THE NODES, DQN(I), I=0,N +********************************************************************* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION QN(0:*), CD(0:*), DQN(0:*) + IF(N .EQ. 0) RETURN + + CALL FCCHGL(N,QN,DQN) + CD(N) = 0.D0 + CD(0) = DQN(1) + IF(N .EQ. 1) GOTO 2 + + DN = DFLOAT(N) + CD(N-1) = 2.D0*DN*DQN(N) + DO 1 K=0,N-2 + KR = N-K-2 + IF(KR .NE. 0) THEN + DK = 2.D0*(DFLOAT(KR)+1.D0) + CD(KR) = CD(KR+2)+DK*DQN(KR+1) + ELSE + CD(0) = .5D0*CD(2)+DQN(1) + ENDIF +1 CONTINUE + +2 CALL FVCHGL(N,CD,DQN) + + RETURN + END +C + SUBROUTINE EDLEGL(N,ET,VN,FU,S1,S2,PA,SO,DSO,D2SO) +************************************************************************ +* COMPUTES THE APPROXIMATE SOLUTION OF A DIRICHLET PROBLEM BY +* COLLOCATION AT THE LEGENDRE GAUSS-LOBATTO NODES +* N = DEGREE OF THE APPROXIMATING POLYNOMIAL +* ET = VECTOR OF THE NODES, ET(I), I=0,N +* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N +* FU = VALUES OF THE RIGHT-HAND SIDE AT THE NODES, FU(I), I=0,N +* S1 = DIRICHLET DATUM AT X=-1 +* S2 = DIRICHLET DATUM AT X=1 +* PA = PARAMETER >0 +* SO = VALUES OF THE APPROXIMATED SOLUTION AT THE NODES, SO(I), I=0,N +* DSO = DERIVATIVE OF THE SOLUTION AT THE NODES, DSO(I), I=0,N +* D2SO= SECOND DERIVATIVE OF THE SOLUTION AT THE NODES, D2SO(I), I=0,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ET(0:*), VN(0:*), FU(0:*), SO(0:*), DSO(0:*), D2SO(0:*) + IF (N .EQ. 0) RETURN + + SO(0) = S1 + SO(N) = S2 + C1 = .5D0*(S2-S1) + DSO(0) = C1 + DSO(N) = C1 + D2SO(0) = 0.D0 + D2SO(N) = 0.D0 + IF (N .EQ. 1) RETURN + + DO 1 J=1,N-1 + SO(J) = .5D0*(S1+S2)+C1*ET(J) +1 CONTINUE + + TH = .58D0 + IM = 42+DLOG(1.D0+PA) + DO 2 IT=1,IM + CALL DELEGL(N,ET,VN,SO,DSO) + CALL DELEGL(N,ET,VN,DSO,D2SO) + DO 3 I=1,N-1 + D2SO(I) = D2SO(I)+FU(I)-PA*SO(I) +3 CONTINUE + + D2SO(N) = 0.D0 + C2 = 2.D0/(ET(1)-ET(0)) + BE = C2/(ET(2)-ET(0)) + DSO(1) = D2SO(1) + C3 = 0.D0 + IF(N .EQ. 2) GOTO 5 + + DO 4 I=1,N-2 + D1 = ET(I)-ET(I-1) + D2 = ET(I+1)-ET(I) + D3 = ET(I+2)-ET(I+1) + C4 = -BE*C3+2.D0/(D1*D2)+PA + D2SO(I) = C4 + BE = 2.D0/(C4*D2*(D2+D3)) + DSO(I+1) = D2SO(I+1)+BE*DSO(I) + C3 = 2.D0/(D2*(D1+D2)) +4 CONTINUE + +5 D2SO(N-1) = -BE*C3+C2/(ET(2)-ET(1))+PA + DSO(N) = 0.D0 + DO 6 I=1,N-1 + IR = N-I + D1 = ET(IR)-ET(IR-1) + D2 = ET(IR+1)-ET(IR) + DSO(IR) = (DSO(IR)+2.D0*DSO(IR+1)/(D2*(D1+D2)))/D2SO(IR) +6 CONTINUE + + DO 7 I=1,N-1 + SO(I) = SO(I)+TH*DSO(I) +7 CONTINUE + +2 CONTINUE + + CALL DELEGL(N,ET,VN,SO,DSO) + CALL DELEGL(N,ET,VN,DSO,D2SO) + + RETURN + END +C + SUBROUTINE DMLEGA(N,NM,CS,DZ,DGA) +************************************************************************ +* COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE +* LEGENDRE GAUSS NODES +* N = DIMENSION OF THE MATRIX +* NM = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT +* CS = VECTOR OF THE ZEROES, CS(I), I=1,N +* DZ = LEGENDRE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N +* DGA = DERIVATIVE MATRIX, DGA(I,J), I=1,N J=1,N +************************************************************************ + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION CS(1), DZ(1), DGA(NM,NM) + DGA(1,1) = 0.D0 + IF (N .LE. 1) RETURN + + DO 1 I=1,N + VI = DZ(I) + ZI = CS(I) + DO 2 J=1,N + IF (I .NE. J) THEN + VJ = DZ(J) + ZJ = CS(J) + DGA(I,J) = VI/(VJ*(ZI-ZJ)) + ELSE + DGA(I,I) = ZI/(1.D0-ZI*ZI) + ENDIF +2 CONTINUE +1 CONTINUE + + RETURN + END