import Test.HUnit
import Text.Printf
import Numeric.Tools.Integration
import Numeric.Tools.Differentiation
import Numeric.Tools.Equation
import Numeric.ApproxEq
----------------------------------------------------------------
-- Integration
----------------------------------------------------------------
-- Functions together with indefinite integral and list of ranges
type FunctionInt = ( Double -> Double -- Function
, Double -> Double -- Integral
, [(Double,Double)] -- List of ranges
, String -- Name
)
-- Test integrator
integratorTest :: String
-> (QuadParam -> (Double,Double) -> (Double -> Double) -> QuadRes)
-> QuadParam
-> FunctionInt
-> [Test]
integratorTest name quad param (f,f',ranges,fname) =
[ case quadRes $ quad param (a,b) f of
Nothing -> TestCase $ assertFailure (printf "%s: convergence for %s failed (%f,%f)" name fname a b)
Just appr ->
let exact = f' b - f' a
in TestCase $ assertBool (printf "%s: poor convergence for %s (%f,%f) %g instead of %g"
name fname a b appr exact)
(eqRelative (quadPrecision param) exact appr)
| (a,b) <- ranges
]
testIntegration :: [Test]
testIntegration = concat
[ integratorTest "Trapeze" quadTrapezoid defQuad =<< [funBlamg,funExp,funLog]
, integratorTest "Simpson" quadSimpson defQuad =<< [funBlamg,funExp,funLog]
, integratorTest "Romberg" quadRomberg defQuad =<< [funBlamg,funExp,funLog]
]
where
funBlamg = ( \x -> x^4 * log(x + sqrt (x*x + 1))
, \x -> 1/5*x^5 * log(x + sqrt (x*x + 1)) - 1/75*sqrt(x*x+1)*(3*x^4 - 4*x^2 + 8)
, [(0,2), (1,3), (-2,3)]
, "x^4·log(x + sqrt(x^2 + 1))"
)
funExp = ( exp, exp
, [(0,2), (1,3), (-2,3)]
, "exp"
)
funLog = ( log
, \x -> x * (log x - 1)
, [(1,2), (0.3,3)]
, "log"
)
----------------------------------------------------------------
--
----------------------------------------------------------------
type FunctionDiff = ( Double -> Double -- Function
, Double -> Double -- Derivative
, [(Double,Double)] -- Points and delta to evaluate
, String -- Name
)
differentiationTest :: String
-> ((Double -> Double) -> Double -> Double -> DiffRes)
-> FunctionDiff
-> [Test]
differentiationTest name diff (f,f',xs,fname) =
[ let DiffRes appr err = diff f h x
exact = f' x
in TestCase $
assertBool
(printf "%s: poor precision for %s, got %g instead of %g" name fname appr exact)
(eqRelative 1e-13 appr exact)
| (x,h) <- xs
]
testDifferentiation :: [Test]
testDifferentiation = concat
[ differentiationTest "richardson" diffRichardson =<< [funSqr,funExp,funLog]
]
where
funSqr = ( \x -> x*x
, \x -> 2*x
, zip ([-10 .. -1]++[1..10]) (repeat 1)
, "square"
)
funExp = ( exp, exp
, zip [-10..10] (repeat 1)
, "exp"
)
funLog = ( log, recip
, map (\x -> (x,x/3)) [0.1,0.2 .. 2]
, "log"
)
testEquation :: [Test]
testEquation =
[ TestCase $ assertBool "Bisection" $ ok (pi/2) (solveBisection 0 (1,2) cos)
, TestCase $ assertBool "Ridders" $ ok (pi/2) $ solveRidders 0 (1,2) cos
, TestCase $ assertBool "Newton" $ ok (pi/2) $ solveNewton 0 (1,2) cos sin
]
where
ok exact (Root x) = within 1 exact x
ok _ _ = False
main :: IO ()
main = do
res <- runTestTT $ TestList $ concat [ testDifferentiation
, testIntegration
, testEquation
]
print res