simpsons_rule Function

public function simpsons_rule(f, dx) result(I)

Performs 1D integration using Simpsons rule

###References * Simpson's rule

Arguments

Type IntentOptional AttributesName
real(kind=Float64), intent(in), dimension(:):: f

Array of equally spaced values

real(kind=Float64), intent(in) :: dx

Spacing between x values

Return Value real(kind=Float64)


Called by

proc~~simpsons_rule~~CalledByGraph proc~simpsons_rule simpsons_rule proc~bt_maxwellian_eb bt_maxwellian_eb proc~bt_maxwellian_eb->proc~simpsons_rule proc~bt_maxwellian_q_n bt_maxwellian_q_n proc~bt_maxwellian_q_n->proc~simpsons_rule proc~bt_maxwellian_n bt_maxwellian_n proc~bt_maxwellian_n->proc~simpsons_rule proc~bt_maxwellian_n_m bt_maxwellian_n_m proc~bt_maxwellian_n_m->proc~simpsons_rule proc~bt_maxwellian_q_n_m bt_maxwellian_q_n_m proc~bt_maxwellian_q_n_m->proc~simpsons_rule interface~bt_maxwellian bt_maxwellian interface~bt_maxwellian->proc~bt_maxwellian_eb interface~bt_maxwellian->proc~bt_maxwellian_q_n interface~bt_maxwellian->proc~bt_maxwellian_n interface~bt_maxwellian->proc~bt_maxwellian_n_m interface~bt_maxwellian->proc~bt_maxwellian_q_n_m

Contents

Source Code


Source Code

function simpsons_rule(f, dx) result(I)
    !+ Performs 1D integration using Simpsons rule
    !+
    !+ ###References
    !+* [Simpson's rule](http://mathworld.wolfram.com/SimpsonsRule.html)
    real(Float64), dimension(:), intent(in) :: f
        !+ Array of equally spaced \(f(x)\) values
    real(Float64), intent(in)               :: dx
        !+ Spacing between x values
    real(Float64)                           :: I

    integer :: s, ii

    s = size(f)
    I = 0.d0
    if(mod(s,2).eq.1) then
        write(*,'(a)') "Length of array must be even"
        return
    endif

    I = f(1)
    do ii=2,s-1
        if(mod(ii,2).eq.1) then
            I = I + 4.0*f(ii)
        else
            I = I + 2.0*f(ii)
        endif
    enddo
    I = I + f(s)
    I = (dx/3.0)*I

end function simpsons_rule