[updated 10.Jul.2005]
Librairie ef3d > Fichier segmqd.f |
SUBROUTINE SEGMQD (NDIM,NUMERO,NBPTQ,POIDQ,COORQ,NSAUT,DEXACT)
Auteur : D.Martin (Fevrier 1988)
Derniere modification : D.Martin (08 juillet 2005)
Version 2
Formules de quadrature sur le segment [0.,1.]
(programme pour que les valeurs des coeff. soit corrects en double precision
pour les compilations implicites en double precision)
-- Variables d'entree --
NDIM dimension d'espace
NUMERO numero du schema de quadrature
(par convention = degre de precision pour les formules de Gauss)
= 100+degre de precision pour les formules de Gauss-Lobatto)
-- Variables en sortie --
NBPTQ nombre de points du schema de quadrature
POIDQ tableau des poids du schema de quadrature
COORQ tableau des coordonnees des points de quadrature
DEXACT degre de precision de la formule (pour impression)
IMPLICIT NONE
INTEGER NDIM,NUMERO,NBPTQ,NSAUT,DEXACT
REAL POIDQ(*),COORQ(NDIM,*)
INTEGER ICQ,IPQ,I
REAL AUX
INTEGER PTQMAX,DEGMAX,DIMMAX
PARAMETER (DEGMAX = 99, PTQMAX=(DEGMAX+1)/2, DIMMAX=(PTQMAX+1)/2)
Poids et points de quadrature de Gauss sur [0,1]
------------------------------------------------
DOUBLE PRECISION PGAUS2,PGAUS3(2),POGAUS(DIMMAX)
DOUBLE PRECISION CGAUS2,CGAUS3(2),COGAUS(DIMMAX)
Formule de Gauss a 2 points, exacte pour P3(0,1)
DATA PGAUS2/0.5D0/
DATA CGAUS2/0.211324865405187D0/
Formule de Gauss a 3 points, exacte pour P5(0,1)
DATA PGAUS3/0.2777777777777778D0,0.4444444444444444D0/
DATA CGAUS3/0.1127016653792585D0,0.5D0/
IF (NUMERO.EQ.-1.OR.NUMERO.EQ.201) THEN
Formule 'nodale' du trapeze exacte pour P1(0,1)
(Newton-Cotes closed type)
NBPTQ=2
DEXACT=1
COORQ(1,1)=1.
COORQ(1,2)=0.
POIDQ(1)=.5
POIDQ(2)=.5
ELSEIF (NUMERO.EQ.-2.OR.NUMERO.EQ.203) THEN
Formule 'nodale' de Simpson exacte pour P3(0,1)
(Newton-Cotes closed type)
NBPTQ=3
DEXACT=3
COORQ(1,1)=1.
COORQ(1,2)=0.
COORQ(1,3)=.5
AUX = 1./6.
POIDQ(1)=AUX
POIDQ(2)=AUX
POIDQ(3)=4.*AUX
ELSEIF (NUMERO.EQ.-3) THEN
Formule 'nodale' de Simpson des 3/8eme exacte pour P3(0,1)
(Newton-Cotes closed type)
NBPTQ=4
DEXACT=3
COORQ(1,1)=1.
COORQ(1,2)=0.
COORQ(1,3)=2./3
COORQ(1,4)=1./3
AUX = .125
POIDQ(1)=AUX
POIDQ(2)=AUX
POIDQ(3)=3*AUX
POIDQ(4)=3*AUX
ELSEIF (NUMERO.EQ.-4) THEN
Formule 'nodale' de Bode exacte pour P5(0,1)
(Newton-Cotes closed type)
NBPTQ=5
DEXACT=5
COORQ(1,1) = 1.
COORQ(1,2) = 0.
COORQ(1,3) = 0.75
COORQ(1,4) = 0.5
COORQ(1,5) = 0.25
AUX = 1./90
POIDQ(1) = 7*AUX
POIDQ(2) = 7*AUX
POIDQ(3) = 32*AUX
POIDQ(4) = 12*AUX
POIDQ(5) = 32*AUX
ELSEIF (NUMERO.EQ.1) THEN
Formule du rectangle exacte pour P1(0.,1.)
NBPTQ=1
DEXACT=1
COORQ(1,1) = .5
POIDQ(1) = 1.
ELSEIF (NUMERO.EQ.3) THEN
Formule de Gauss a 2 points exacte pour P3(0.,1.)
[--1-----2--]
NBPTQ=(NUMERO+1)/2
DEXACT=NUMERO
IF (NSAUT.EQ.0) NSAUT=1
COORQ(1,1) = CGAUS2
POIDQ(1) = .5
COORQ(1,1+NSAUT) = 1.-CGAUS2
POIDQ(2) = .5
ELSEIF (NUMERO.EQ.5) THEN
Formule de Gauss a 3 points exacte pour P5(0.,1.)
[-1---2---3-]
NBPTQ=(NUMERO+1)/2
DEXACT=NUMERO
IF (NSAUT.LE.0) NSAUT=1
AUX = 5./18.
COORQ(1,1) = CGAUS3(1)
POIDQ(1) = AUX
COORQ(1,1+NSAUT) = .5
POIDQ(2) = 4./9.
COORQ(1,1+2*NSAUT) = 1.-CGAUS3(1)
POIDQ(3) = AUX
ELSEIF (MOD (NUMERO,2).EQ.1.AND.NUMERO.LE.DEGMAX) THEN
Formule de Gauss a n points exacte pour P{2n-1}(0.,1.)
DEXACT=NUMERO
NBPTQ=(DEXACT+1)/2
CALL GAUSLE (NBPTQ,COGAUS,POGAUS)
IF (NSAUT.LE.0) NSAUT=1
ICQ=1
IPQ=1
DO 111 I=1,NBPTQ/2
COORQ(1,ICQ) = COGAUS(I)
POIDQ(IPQ) = POGAUS(I)
ICQ=ICQ+NSAUT
IPQ=IPQ+1
111 CONTINUE
IF (MOD (NBPTQ,2).EQ.1) THEN
COORQ(1,ICQ) = .5
POIDQ(IPQ) = POGAUS((NBPTQ+1)/2)
ICQ=ICQ+NSAUT
IPQ=IPQ+1
ENDIF
IPQ=NBPTQ+1
ICQ=NBPTQ*NSAUT+1
DO 112 I=1,NBPTQ/2
ICQ=ICQ-NSAUT
IPQ=IPQ-1
COORQ(1,ICQ) = 1.-COGAUS(I)
POIDQ(IPQ) = POGAUS(I)
112 CONTINUE
ELSEIF (MOD (NUMERO,2).EQ.1.AND.NUMERO.GT.DEGMAX+1) THEN
Formule de Gauss-Lobatto a n points exacte pour P_{2n-3}(0.,1.)
DEXACT=NUMERO-DEGMAX-1
NBPTQ=(DEXACT+3)/2
IF (NSAUT.LE.0) NSAUT=1
ICQ=1+NSAUT
IPQ=2
COORQ(1,1) = 1.
POIDQ(1) = 1./(NBPTQ*(NBPTQ-1))
COORQ(1,ICQ) = 0.
POIDQ(IPQ) = POIDQ(1)
IF (NBPTQ.LE.2) RETURN
CALL GAUSLO (NBPTQ,COGAUS,POGAUS)
DO 122 I=1,(NBPTQ-2)/2
ICQ=ICQ+NSAUT
IPQ=IPQ+1
COORQ(1,ICQ) = 1.-COGAUS(I)
POIDQ(IPQ) = POGAUS(I)
122 CONTINUE
IF (MOD (NBPTQ,2).EQ.1) THEN
ICQ=ICQ+NSAUT
IPQ=IPQ+1
COORQ(1,ICQ) = .5
POIDQ(IPQ) = POGAUS((NBPTQ-1)/2)
ENDIF
DO 121 I=(NBPTQ-2)/2,1,-1
ICQ=ICQ+NSAUT
IPQ=IPQ+1
COORQ(1,ICQ) = COGAUS(I)
POIDQ(IPQ) = POGAUS(I)
121 CONTINUE
ELSE
CALL SEXXQD (NDIM,NUMERO,NBPTQ,POIDQ,COORQ,DEXACT)
ENDIF
END !Segmqd
segmqd est appelé dans (4 procédures)