XProfan

WindowTitle "QANC8 Integration"
'{----------------------------------------------------------
' Programm zur schnellen näherungsweisen Berechnung von Integralen
' mittels der https://de.wikipedia.org/wiki/Newton-Cotes-Formeln .
' Im Artikel wird die Formelordnung zwischen 2 und 6 erläutert.
' Das nachstehende Programm verwendet eine Formel 8. Ordnung mit
' vom Anwender vorgebbarer absoluter und relativer Genauigkeit.
' Von Fortran-90 nach XProfan11 2016-10 by P.Specht, Vienna/Austria/EU
' Keine wie auch immer geartete Gewähr! Möglicherweise Rechte Dritter!
' ------------------------------------------------------------------
' Program to integrate a user-defined function f(x) from x1 to x2 by
' the QANC8 subroutine with control of abs. and relative precisions
' ------------------------------------------------------------------
' SAMPLE RUN:
'  (Integrate function cos(x) - 2 sin(x) from x=0 to x=1,
'   with a precision <= 1e-10):
'
'  Integral Value = -7.792440345582408E-002
'  Estimated error =  1.041851518709824E-017
'  Error code =  0.000000000000000E+000
'  Number of function evaluations =          33
' ------------------------------------------------------------------
' Reference: From Numath Library By Tuan Dang Trong in Fortran 77
'                              F90 Release 1.0 By J-P Moreau, Paris
' F90 translated to XProfan 11.2a by P.Specht, Vienna /Austria.
' As based on piecewise polynomial approximation, quanc8 is not
' designed to handle certain kinds of integrals (e.g. functions f(x)
' where derivatives < 10th order are unbound or do not exist).
'}---------------------------------------------------------------
' PROGRAM TQANC8  'Test Quanc8
Windowstyle 24:Cls rgb(200,200,220):Font 2
Declare  AERROR!,CODE![0],Error![0],RERROR!,X1!,X2!,Valu![0],NBRF&[0]
X1!=0
X2!=1
AERROR!=Val("1e-9")
RERROR!=Val("1e-10")
QANC8(X1!,X2!,AERROR!,RERROR!,Valu![],Error![],NBRF&[],CODE![])
print "\n Integral  Value =", format$("%e",Valu![0])
Print "\n Estimated error =", format$("%e",Error![0])
Print "\n Error Code =", format$("%g",CODE![0])
Print "\n Number of function evaluations =", format$("%u",NBRF&[0])
Beep
WaitInput
End
'--------------------------------------------------------------
'  Zu integrierende Funktion:
'--------------------------------------------------------------

Proc FCT : Parameters X!

    ' Für ein allgemeintes Intervall [x1,x2] sind
    ' als Stützstellen x = x1 +(x2-x1)*x zu nehmen!
    Declare fct!
    FCT!=Cos(X!)-2*Sin(X!)
    Return FCT!

EndProc

proc DMAX1 :parameters z1!,z2! :return if(z1!>z2!,z1!,z2!)

endproc

'{--------------------------------------------------------------
' QUICK APPROXIMATION using NEWTON COTES of Order 8:
' Proc  QANC8 :Parameters A!,B!,AERR!,RERR!,RES!,ERR!,NBF&,FLG!)
'
'     INTEGRATE A REAL FUNCTION FCT(X) FROM X=A TO X=B, WITH
'     GIVEN ABSOLUTE AND RELATIVE PRECISIONS, AERR, RERR.
'     INPUTS:
'     FCT     EXTERNAL USER-DEFINED FUNCTION FOR ANY X VALUE
'             IN INTERVAL (A,B)
'     A,B     LIMITS OF INTERVAL
'     AERR,RERR   RESPECTIVELY ABSOLUTE ERROR AND RELATIVE ERROR
'                 REQUIRED BY USER
'     OUTPUTS:
'     RES     VALUE OF INTEGRAL
'     ERR     ESTIMATED ERROR
'     NBF     NUMBER OF NECESSARY FCT(X) EVALUATIONS
'     FLG     INDICATOR
'             = 0.0       CORRECT RESULT
'             = NNN.RRR   NO CONVERGENCE DU TO A SINGULARITY.
'             THE SINGULAR POINT ABCISSA IS GIVEN BY FORMULA:
'             XS = B-.RRR*(B-A)
' Ref.: FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL
'       COMPUTATIONS. PRENTICE-HALL, INC.
' ------------------------------------------------------------
'}    IMPLICIT REAL *8 (A-H,O-Z)

Proc  QANC8 :Parameters A!,B!,AERR!,RERR!,RES![],ERR![],NBF&[],FLG![]

    Declare LMIN&,LMAX&,LOUT&,NMAX&,NFIN&,W0!,W1!,W2!,W3!,W4!
    Declare Cor!,Sum!,L&,NIM&,X0!,QP!,PAS!,PAS1!,I&,J&,F0!
    Declare QL!,QN!,QD!,ERR1!,Tol1!,TEMP!
    Declare QR![31],F![16],X![16],FS![8,30],XS![8,30]
    LMIN& = 1
    LMAX& = 30
    LOUT& = 6
    NMAX& = 5000
    NFIN& = NMAX&-8*(LMAX&-LOUT&+2^(LOUT&+1))
    W0!  =   3956/14175
    W1!  =  23552/14175
    W2!  =  -3712/14175
    W3!  =  41984/14175
    W4!  = -18160/14175
    FLG![0] = 0
    RES![0] = 0
    COR! = 0
    ERR![0] = 0
    SUM! = 0
    NBF&[0] = 0
    Case  A!=B!: Return
    L& = 0
    NIM& = 1
    X0!  = A!
    X![16] = B!
    QP! = 0
    F0!   = FCT(X0!)
    PAS1!  = (B!-A!)/16
    X![8]  = (X0!+X![16])/2
    X![4]  = (X0!+X![8])/2
    X![12] = (X![8]+X![16])/2
    X![2]  = (X0!+X![4])/2
    X![6]  = (X![4]+X![8])/2
    X![10] = (X![8]+X![12])/2
    X![14] = (X![12]+X![16])/2

    Whileloop 2,16,2:j&=&Loop

        F![J&] = FCT(X![J&])

    EndWhile

    NBF&[0] = 9
    L30:
    X![1]  = (X0!+X![2])/2
    F![1] = FCT(X![1])

    WhileLoop 3,15,2:J&=&Loop

        X![J&]  = (X![J&-1]+X![J&+1])/2
        F![J&] = FCT(X![J&])

    EndWhile

    L35:
    NBF&[0] = NBF&[0]+8
    PAS! = (X![16]-X0!)/16
    QL!  = (W0!*(F0!+F![8])+W1!*(F![1]+F![7])+\
    W2!*(F![2]+F![6])+W3!*(F![3]+F![5])+W4!*F![4])*PAS!
    QR![L&+1] = (W0!*(F![8]+F![16])+W1!*(F![9]+F![15])+\
    W2!*(F![10]+F![14])+W3!*(F![11]+F![13])+W4!*F![12])*PAS!
    QN! = QL! + QR![L&+1]
    QD! = QN! - QP!
    SUM! = SUM! + QD!
    ERR1! = Abs(QD!)/1023
    TOL1! = DMAX1(AERR!, RERR!*Abs(SUM!)) * (PAS!/PAS1!)
    Case L&<LMIN&:Goto "L50"
    Case L&>=LMAX&:Goto "L62"
    Case NBF&[0]>NFIN&:Goto "L60"
    Case ERR1!<=TOL1!:Goto "L70"
    L50:
    NIM& = 2*NIM&
    L& = L&+1

    WhileLoop 8:i&=&Loop

        FS![I&,L&] = F![I&+8]
        XS![I&,L&] = X![I&+8]

    EndWhile

    L52:
    QP! = QL!

    WhileLoop 8:i&=&Loop

        F![18-2*I&] = F![9-I&]
        X![18-2*I&] = X![9-I&]

    EndWhile

    L55:
    Goto "L30"
    L60:
    NFIN& = 2*NFIN&
    LMAX& = LOUT&
    FLG![0] = FLG![0] + (B!-X0!)/(B!-A!)
    Goto "L70"
    L62:
    FLG![0] = FLG![0] + 1
    L70:
    RES![0] = RES![0] + QN!
    ERR![0] = ERR![0] + ERR1!
    COR! = COR! + QD!/1023
    L72:
    ' (A!-H!,O!-Z!)
    Case NIM&=(int(NIM&/2)*2):Goto "L75"'Gerade
    NIM& = NIM&/2
    L& = L&-1
    Goto "L72"
    L75:
    NIM& = NIM&+1
    Case L&<=0:Goto "L80"
    QP! = QR![L&]
    X0! = X![16]
    F0! = F![16]

    WhileLoop 8:i&=&Loop

        F![2*I&] = FS![I&,L&]
        X![2*I&] = XS![I&,L&]

    Endwhile

    Goto "L30"
    L80:
    RES![0] = RES![0] + COR!
    Case ERR![0]=0: RETURN
    L82:
    TEMP! = Abs(RES![0]) + ERR![0]
    Case TEMP!<>Abs(RES![0]): RETURN
    ERR![0] = 2*ERR![0]
    Goto "L82"

EndProc