aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorWojtek Kosior <kwojtus@protonmail.com>2019-06-21 13:08:31 +0200
committerWojtek Kosior <kwojtus@protonmail.com>2019-06-21 13:08:31 +0200
commit7d28fe599e1a5498a7e86dc04fa207cab11559ac (patch)
treecf9695c58ece1f808a9c67f2beff751d7de0720f
parent63c2faecd587b08b4a481a0da8b821edd58db572 (diff)
downloadfortran-assignment3-7d28fe599e1a5498a7e86dc04fa207cab11559ac.tar.gz
fortran-assignment3-7d28fe599e1a5498a7e86dc04fa207cab11559ac.zip
newton-cotes integration of entire interval at once
-rw-r--r--src/quadratures.f9049
1 files changed, 45 insertions, 4 deletions
diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 90bba6e..a3f869f 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -21,13 +21,22 @@ MODULE quadratures
END INTERFACE
INTERFACE
+ FUNCTION quadrature(qbeg, qend, fun) result(val)
+ IMPLICIT none
+ real(kind=8), intent(in) :: qbeg, qend
+ procedure(funint) :: fun
+ real(kind=8) :: val
+ END FUNCTION quadrature
+ END INTERFACE
+
+ INTERFACE
FUNCTION funint(x) result(y)
IMPLICIT none
real(kind=8), intent(in) :: x
real(kind=8) :: y
END FUNCTION funint
END INTERFACE
-
+
CONTAINS
FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
@@ -38,19 +47,51 @@ CONTAINS
integer(kind=4), intent(in) :: p
real(kind=8) :: val
+ procedure(quadrature), pointer :: quad
+
SELECT CASE (p)
CASE (:-1)
STOP "negative interpolationg polynomial order passed"
CASE (0)
- val = 0 ! rectangle(ibeg, iend, fun)
+ quad => rectangle
CASE (1)
- val = 0 ! trapeze(ibeg, iend, fun)
+ quad => trapeze
CASE (2)
- val = 0 ! simpson_1_third(ibeg, iend, fun)
+ quad => simpson_1_third
CASE default
STOP "Newton-Cotes quadratures only implemented for order < 3"
END SELECT
+
+ ! compute integral using quadrature pointed by quad
+ val = quad(ibeg, iend, fun)
END FUNCTION newton_cotes
+ FUNCTION rectangle(qbeg, qend, fun) result(val)
+ IMPLICIT none
+ real(kind=8), intent(in) :: qbeg, qend
+ procedure(funint) :: fun
+ real(kind=8) :: val
+
+ val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
+ END FUNCTION rectangle
+
+ FUNCTION trapeze(qbeg, qend, fun) result(val)
+ IMPLICIT none
+ real(kind=8), intent(in) :: qbeg, qend
+ procedure(funint) :: fun
+ real(kind=8) :: val
+
+ val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
+ END FUNCTION trapeze
+ FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
+ IMPLICIT none
+ real(kind=8), intent(in) :: qbeg, qend
+ procedure(funint) :: fun
+ real(kind=8) :: val
+
+ val = (qend - qbeg) * (1/6.0) * &
+ (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
+ END FUNCTION simpson_1_third
+
END MODULE quadratures