# calculus.test --
#    Test cases for the Calculus package
#
# This file contains a collection of tests for one or more of the Tcllib
# procedures.  Sourcing this file into Tcl runs the tests and
# generates output for errors.  No output means no errors were found.
#
# Copyright (c) 2002, 2003, 2004 by Arjen Markus.
# Copyright (c) 2004 by Kevin B. Kenny
# All rights reserved.
#
# RCS: @(#) $Id: calculus.test,v 1.17 2008/06/26 20:04:23 arjenmarkus Exp $

# -------------------------------------------------------------------------

source [file join \
	[file dirname [file dirname [file join [pwd] [info script]]]] \
	devtools testutilities.tcl]

testsNeedTcl     8.4
testsNeedTcltest 2.1

support {
    useLocal math.tcl        math
    useLocal interpolate.tcl math::interpolate
}
testing {
    useLocal calculus.tcl math::calculus
}

# -------------------------------------------------------------------------

package require log
log::lvSuppress notice

# -------------------------------------------------------------------------

namespace eval ::math::calculus::test {

namespace import ::tcltest::test
namespace import ::math::calculus::*

set prec $::tcl_precision
set ::tcl_precision 17

proc matchNumbers {expected actual} {
    set match 1
    foreach a $actual e $expected {
        if {abs($a-$e) > 0.1e-4} {
            set match 0
            break
        }
    }
    return $match
}

customMatch numbers matchNumbers

#
# Simple test functions - exact result predictable!
#
proc const_func { x } {
   return 1
}
proc linear_func { x } {
   return $x
}
proc downward_linear { x } {
   return [expr {100.0-$x}]
}

#
# Test the Integral proc
#
test "Integral-1.0" "Integral of constant function" {
   integral 0 100 100 const_func
} 100.0

test "Integral-1.1" "Integral of linear function" {
   integral 0 100 100 linear_func
} 5000.0

test "Integral-1.2" "Integral of downward linear function" {
   integral 0 100 100 downward_linear
} 5000.0

test "Integral-1.3" "Integral of expression" {
   integralExpr 0 100 100 {100.0-$x}
} 5000.0


proc const_func2d { x y } {
   return 1
}
proc linear_func2d { x y } {
   return $x
}

test "Integral2D-1.0" "Integral of constant 2D function" {
   integral2D { 0 100 10 } { 0 50 1 } const_func2d
} 5000.0
test "Integral2D-1.1" "Integral of constant 2D function (different step)" {
   integral2D { 0 100 1 } { 0 50 1 } const_func2d
} 5000.0
test "Integral2D-1.2" "Integral of linear 2D function" {
   integral2D { 0 100 10 } { 0 50 1 } linear_func2d
} 250000.0


proc const_func3d { x y z } {
   return 1
}
proc linear_func3d { x y z } {
   return $x
}

test "Integral3D-1.0" "Integral of constant 2D function" {
   integral3D { 0 100 10 } { 0 50 1 } { 0 50 1 } const_func3d
} 250000.0
test "Integral3D-1.1" "Integral of constant 2D function (different step)" {
   integral3D { 0 100 1 } { 0 50 1 } { 0 50 1 } const_func3d
} 250000.0
test "Integral3D-1.2" "Integral of linear 2D function" {
   integral3D { 0 100 10 } { 0 50 1 } { 0 50 1 } linear_func3d
} 12500000.0

proc f2d_1 {x y} {
    return 1
}
proc f2d_x {x y} {
    return $x
}
proc f2d_y {x y} {
    return $y
}
proc f2d_x2 {x y} {
    return [expr {$x*$x}]
}
proc f2d_y2 {x y} {
    return [expr {$y*$y}]
}

test "Integral2D-2.0" "Integrals of 2D functions - accurate" -match numbers -body {
    set result {}
    foreach f {f2d_1 f2d_x f2d_y f2d_x2 f2d_y2} {
        lappend  result [::math::calculus::integral2D_accurate {-1 1 1} {-1 1 1} $f]
    }
    return $result
} -result {4.0 0.0 0.0 1.333333333 1.333333333}


proc f3d_1 {x y z} {
    return 1
}
proc f3d_x {x y z} {
    return $x
}
proc f3d_y {x y z} {
    return $y
}
proc f3d_z {x y z} {
    return $z
}
proc f3d_x2 {x y z} {
    return [expr {$x*$x}]
}
proc f3d_y2 {x y z} {
    return [expr {$y*$y}]
}
proc f3d_z2 {x y z} {
    return [expr {$z*$z}]
}

test "Integral2D-2.0" "Integrals of 2D functions - accurate" -match numbers -body {
    set result {}
    foreach f {f3d_1 f3d_x f3d_y f3d_z f3d_x2 f3d_y2 f3d_z2} {
        lappend  result [::math::calculus::integral3D_accurate {-1 1 1} {-1 1 1} {-1 1 1} $f]
    }
    return $result
} -result {8.0 0.0 0.0 0.0 2.666666667 2.666666667 2.666666667}


#
# Test cases: yet to be brought into the tcltest form!
#

# xvec should one long!
proc const_func { t xvec } { return 1.0 }

# xvec should be two long!
proc dampened_oscillator { t xvec } {
   set x  [lindex $xvec 0]
   set x1 [lindex $xvec 1]
   return [list $x1 [expr {-$x1-$x}]]
}

foreach method {eulerStep heunStep rungeKuttaStep} {
   log::log notice "Method: $method"

   set xvec   0.0
   set t      0.0
   set tstep  1.0
   for { set i 0 } { $i < 10 } { incr i } {
      set result [$method $t $tstep $xvec const_func]
      log::log notice "Result ($t): $result"
      set t      [expr {$t+$tstep}]
      set xvec   $result
   }

   set xvec   { 1.0 0.0 }
   set t      0.0
   set tstep  0.1
   for { set i 0 } { $i < 20 } { incr i } {
      set result [$method $t $tstep $xvec dampened_oscillator]
      log::log notice "Result ($t): $result"
      set t      [expr {$t+$tstep}]
      set xvec   $result
   }
}

#
# Boundary value problems:
#
proc coeffs { x } { return {1.0 0.0 0.0} }
proc forces { x } { return 0.0 }

log::log notice [boundaryValueSecondOrder coeffs forces {0.0 1.0} {100.0 0.0} 10]
log::log notice [boundaryValueSecondOrder coeffs forces {0.0 0.0} {100.0 1.0} 10]

#
# Determining the root of an equation
# use simple functions
#
proc func  { x } { expr {$x*$x-1.0} }
proc deriv { x } { expr {2.0*$x} }

test "NewtonRaphson-1.0" "Result should be 1" {
   set result [newtonRaphson func deriv 2.0]
   if { abs($result-1.0) < 0.0001 } {
      set answer 1
   }
} 1
test "NewtonRaphson-1.1" "Result should be -1" {
   set result [newtonRaphson func deriv -0.5]
   if { abs($result+1.0) < 0.0001 } {
      set answer 1
   }
} 1

proc func2  { x } { expr {$x*exp($x)-1.0} }
proc deriv2 { x } { expr {exp($x)+$x*exp($x)} }

test "NewtonRaphson-2.1" "Result should be nearly 0.56714" {
   set result [newtonRaphson func2 deriv2 2.0]
   if { abs($result-0.56714) < 0.0001 } {
      set answer 1
   }
} 1

test "NewtonRaphson-2.2" "Result should be nearly 0.56714" {
   set result [newtonRaphson func2 deriv2 -0.5]
   if { abs($result-0.56714) < 0.0001 } {
      set answer 1
   }
} 1

proc checkout { expr integrator a b target } {
    set problems {}
    proc g x [list expr $expr]
    set cmd $integrator
    lappend cmd g $a $b
    foreach { s error } [eval $cmd] break
    set diff [expr { abs( $s - $target ) }]
    if { $diff > 1.0e-6 * $target && $diff > 1.0e-10 } {
	append problems \n  "error underestimated!" \
	    \n "f =" $expr ", a=" $a ", b=" $b \
	    \n "machinery = " $integrator "," \
	    \n "estimated " $error " actual " $diff
    }
    return $problems
}

test romberg-1.1 {simple integral} {
    checkout { pow( $x, 16 ) } romberg -1. 1. [expr { 2. / 17. }]
} {}
test romberg-1.2 {simple integral} {
    checkout { exp( -$x * $x / 2. ) / sqrt( 2. * 3.1415926535897932 ) } \
	romberg -1. 1. 0.68268949213708590
} {}
test romberg-1.3 {simple integral} {
    checkout { sin($x) } romberg 0 3.1415926535897932 2.0
} {}

test romberg-1.4  { Singularity where limit exists } {
    checkout { sin($x)/$x } romberg 0 3.1415926535897932 1.8519370519824662
} {}

test romberg-1.5 { Parameter error } {
    catch {romberg irrelevant 0 1 -degree} result
    set result
} "wrong \# args, should be \"romberg f x1 x2 ?-option value?...\""

test romberg-1.6 { Parameter error } {
    catch {romberg irrelevant 0 1 -bad flag} result
    set result
} "unknown option \"-bad\", should be -abserror, -degree, -relerror, or\
   -maxiter"

test romberg-1.7 { Max iterations exceeded } \
    -setup {
	proc f x { expr { pow($x,4) } }
    } \
    -body {
	foreach { value error } [romberg f -1. 1. -degree 1 -maxiter 3 ] break
	expr { abs($value - 0.4) < $error }
    } \
    -cleanup {
	rename f {}
    } \
    -result 1

test romberg-1.8 {Bad param} {
    catch {romberg irrelevant 0 1 -degree bad} result
    set result
} {expected an integer but found "bad"}

test romberg-1.9 {Bad param} {
    catch {romberg irrelevant 0 1 -degree 0} result
    set result
} {-degree must be positive}

test romberg-1.10 {Bad param} {
    catch {romberg irrelevant 0 1 -maxiter bad} result
    set result
} {expected an integer but found "bad"}

test romberg-1.11 {Bad param} {
    catch {romberg irrelevant 0 1 -maxiter 0} result
    set result
} {-maxiter must be positive}

test romberg-1.12 {Bad param} {
    catch {romberg irrelevant 0 1 -abserror bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-1.13 {Bad param} {
    catch {romberg irrelevant 0 1 -abserror 0.} result
    set result
} {-abserror must be positive}

test romberg-1.14 {Bad param} {
    catch {romberg irrelevant 0 1 -relerror bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-1.15 {Bad param} {
    catch {romberg irrelevant 0 1 -relerror 0.} result
    set result
} {-relerror must be positive}

test romberg-1.16 {Bad limit } {
    catch {romberg irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-1.17 {Bad limit} {
    catch {romberg irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-2.1 {Integral over half-infinite interval} {
    checkout { exp( -$x * $x / 2. ) / sqrt( 2. * 3.1415926535897932 ) } \
	romberg_infinity -30. -1. 0.15865525393145705
} {}
test romberg-2.2 {Integral over half-infinite interval} {
    checkout { exp( -$x * $x / 2. ) / sqrt( 2. * 3.1415926535897932 ) } \
	romberg_infinity 1. 30. 0.15865525393145705
} {}
test romberg-2.3 {Integral over half-infinite interval} {
    checkout { exp( $x )  } romberg_infinity -1.e38 -1. [expr { exp(-1.) }]
} {}
test romberg-2.4 {Parameter error} {
    catch {romberg_infinity irrelevant -1.e38 2.} result
    set result
} {limits of integration have opposite sign}

test romberg-2.5 {Bad limit } {
    catch {romberg_infinity irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-2.6 {Bad limit} {
    catch {romberg_infinity irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-3.1 {Square root singularity at the upper bound} {
    checkout { sqrt( 1.0 / ( 1.0 - $x ) ) } romberg_sqrtSingUpper 0. 1. 2.
} {}

test romberg-3.2 \
    {Square root singularity in the derivative at the upper bound} {
	checkout { 4. * sqrt( 1.0 - $x * $x ) } romberg_sqrtSingUpper 0. 1. \
	    3.1415926535897932
    } {}

test romberg-3.3 {Square root singularity at the lower bound} {
    checkout { 1.0 / sqrt($x) } romberg_sqrtSingLower 0. 4. 4.
} {}

test romberg-3.4 \
    {Square root singularity in the derivative at the lower bound} {
	checkout { 4. * sqrt( 1.0 - $x * $x ) } romberg_sqrtSingLower -1. 0. \
	    3.1415926535897932
    } {}

test romberg-3.5 {Bad limit } {
    catch {romberg_sqrtSingUpper irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-3.6 {Bad limit} {
    catch {romberg_sqrtSingUpper irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-3.7 {Bad limits} {
    catch {romberg_sqrtSingUpper irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-3.8 {Bad limit } {
    catch {romberg_sqrtSingLower irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-3.9 {Bad limit} {
    catch {romberg_sqrtSingLower irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-3.10 {Bad limits} {
    catch {romberg_sqrtSingLower irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-4.1 {Power law singularity at the lower bound} {
    checkout { 1.0 / sqrt($x) } [list romberg_powerLawLower 0.5] 0. 4. 4.
} {}

test romberg-4.2 \
    {Power law signularity in the derivative at the lower bound.} {
	checkout { sqrt( sqrt( $x ) ) } \
	    [list romberg_powerLawLower 0.75] 0. 1. 0.8
    } {}

test romberg-4.3 {Power law singularity at the upper bound} {
    checkout { 1.0 / sqrt(4.0 - $x) } \
	[list romberg_powerLawUpper 0.5] 0. 4. 4.
} {}

test romberg-4.4 \
    {Power law singularity in the derivative at the upper bound} {
	checkout { sqrt( sqrt( -$x ) ) } \
	    [list romberg_powerLawUpper 0.75] -1. 0. 0.8
    } {}

test romberg-4.5 {Bad limit } {
    catch {romberg_powerLawUpper 0.5 irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-4.6 {Bad limit} {
    catch {romberg_powerLawUpper 0.5 irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-4.7 {Bad limits} {
    catch {romberg_powerLawUpper 0.5 irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-4.8 {Bad limit } {
    catch {romberg_powerLawLower 0.5 irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-4.9 {Bad limit} {
    catch {romberg_powerLawLower 0.5 irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-4.10 {Bad limits} {
    catch {romberg_powerLawLower 0.5 irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-4.11 {Bad gamma} {
    catch {romberg_powerLawUpper bad irrelevant 1 0} result
    set result
} {expected a floating-point number but found "bad"}
test romberg-4.12 {Bad gamma} {
    catch {romberg_powerLawUpper 0. irrelevant 1. 0.} result
    set result
} {gamma must lie in the interval (0,1)}
test romberg-4.13 {Bad gamma} {
    catch {romberg_powerLawUpper 1. irrelevant 1. 0.} result
    set result
} {gamma must lie in the interval (0,1)}
test romberg-4.14 {Bad gamma} {
    catch {romberg_powerLawLower bad irrelevant 1 0} result
    set result
} {expected a floating-point number but found "bad"}
test romberg-4.15 {Bad gamma} {
    catch {romberg_powerLawLower 0. irrelevant 1. 0.} result
    set result
} {gamma must lie in the interval (0,1)}
test romberg-4.16 {Bad gamma} {
    catch {romberg_powerLawLower 1. irrelevant 1. 0.} result
    set result
} {gamma must lie in the interval (0,1)}

test romberg-5.1 {Function that decays exponentially} {
    checkout { exp( -$x * $x / 2. ) / sqrt( 2. * 3.1415926535897932 ) } \
	romberg_expUpper 1. 100. 0.15865525393145705
} {}

test romberg-5.2 {Function that grows exponentially} {
    checkout { exp( -$x * $x / 2. ) / sqrt( 2. * 3.1415926535897932 ) } \
	romberg_expLower -100. -1. 0.15865525393145705
} {}

test romberg-5.3 {Bad limit } {
    catch {romberg_sqrtSingUpper irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-5.4 {Bad limit} {
    catch {romberg_sqrtSingUpper irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-5.5 {Bad limits} {
    catch {romberg_sqrtSingUpper irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-5.6 {Bad limit } {
    catch {romberg_sqrtSingLower irrelevant bad 1} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-5.7 {Bad limit} {
    catch {romberg_sqrtSingLower irrelevant 0 bad} result
    set result
} {expected a floating-point number but found "bad"}

test romberg-5.8 {Bad limits} {
    catch {romberg_sqrtSingLower irrelevant 1 0} result
    set result
} {limits of integration out of order}

test romberg-6.1 {Fancy integration} \
    -setup {
	proc v {f u} {
	    set x [expr { sin($u) }]
	    set cmd $f; lappend cmd $x; set y [eval $cmd]
	    return [expr { $y * cos($u) }]
	}
	proc romberg_sine { f a b args } {
	    set f [lreplace $f 0 0 \
		       [uplevel 1 [list namespace which [lindex $f 0]]]]
	    set f [list v $f]
	    return [eval [linsert $args 0 \
			      romberg $f \
			      [expr { asin($a) }] [expr { asin($b) }]]]
	}
    } \
    -body {
	checkout { exp($x) / sqrt( 1. - $x * $x ) } romberg_sine -1. 1. \
	    3.97746326
    } \
    -cleanup {
	rename v {}
	rename romberg_sine {}
    } \
    -result {}


proc matchNumbers {expected actual} {
   set match 1
   foreach a $actual e $expected {
      if {abs($a-$e) > 0.1e-6} {
         set match 0
         break
      }
   }
   return $match
}

customMatch numbers matchNumbers

proc ::f1 {x} {expr {1.0-$x}}
proc ::f2 {x} {expr {1.0-$x*$x}}
proc ::f3 {x} {expr {cos($x)}}

test "regula-1.0" "Zero of linear function" \
   -match numbers -body {
   set x1 [::math::calculus::regula_falsi ::f1 0.0 5.0]
} -result 1.0

test "regula-1.1" "Zero of quadratic function" \
   -match numbers -body {
   set x1 [::math::calculus::regula_falsi ::f2 0.0 5.0]
} -result 0.99909822

test "regula-1.2" "Zero of quadratic function (more accurate)" \
   -match numbers -body {
   set x1 [::math::calculus::regula_falsi ::f2 0.0 5.0 1.0e-6]
} -result 0.99999305

test "regula-1.3" "Zero of cosine" \
   -match numbers -body {
   set x1 [::math::calculus::regula_falsi ::f3 0.0 3.0]
} -result 1.5707963

test "regula-2.1" "Negative relative error" \
   -match glob -body {
   set x1 [::math::calculus::regula_falsi ::f1 0.0 3.0 -1.0e-4]
} -result "Relative *" -returnCodes error

test "regula-2.2" "Invalid interval" \
   -match glob -body {
   set x1 [::math::calculus::regula_falsi ::f3 0.0 5.0]
} -result "Interval must be *" -returnCodes error

test "solveTriDiagonal-1.0" "Solve tridiagonal system" \
   -match numbers -body {
   set x [::math::calculus::solveTriDiagonal {3 3} {1 1 1} {2 2} {1 0 0}]
} -result [list [expr {5.0/11.0}] [expr {3.0/11.0}] [expr {-9.0/11.0}]]

# End of test cases
testsuiteCleanup
	
set ::tcl_precision $prec

testsuiteCleanup
}

namespace delete ::math::calculus::test

# Local Variables:
# mode: tcl
# End: