summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2012-02-25 12:40:49 +0100
committerAlberto Ruiz <aruiz@um.es>2012-02-25 12:40:49 +0100
commiteab231907c3d11651e48e15510f4510fd6c77450 (patch)
treeb133814d7b82e230cc11b9eb1e50275aa95e02ae
parent9322df983400894625105f66e36e5718329c0053 (diff)
merge changelog
-rw-r--r--CHANGES.md3
-rw-r--r--examples/ode.hs16
-rw-r--r--lib/Numeric/GSL/ODE.hs47
-rw-r--r--lib/Numeric/GSL/gsl-aux.c62
4 files changed, 87 insertions, 41 deletions
diff --git a/CHANGES.md b/CHANGES.md
index de466e5..10be500 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -3,6 +3,9 @@
3 3
4- integration over infinite intervals 4- integration over infinite intervals
5 5
6- msadams and msbdf methods for ode
7
8
60.13.0.0 90.13.0.0
7-------- 10--------
8 11
diff --git a/examples/ode.hs b/examples/ode.hs
index 082c46c..dc6e0ec 100644
--- a/examples/ode.hs
+++ b/examples/ode.hs
@@ -1,6 +1,9 @@
1{-# LANGUAGE ViewPatterns #-}
1import Numeric.GSL.ODE 2import Numeric.GSL.ODE
2import Numeric.LinearAlgebra 3import Numeric.LinearAlgebra
3import Graphics.Plot 4import Graphics.Plot
5import Debug.Trace(trace)
6debug x = trace (show x) x
4 7
5vanderpol mu = do 8vanderpol mu = do
6 let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)] 9 let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
@@ -32,3 +35,16 @@ main = do
32 harmonic 1 0.1 35 harmonic 1 0.1
33 kepler 0.3 60 36 kepler 0.3 60
34 kepler 0.4 70 37 kepler 0.4 70
38 vanderpol' 2
39
40-- example of odeSolveV with jacobian
41vanderpol' mu = do
42 let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)]
43 jac t (toList->[x,v]) = (2><2) [ 0 , 1
44 , -1-2*x*v*mu, mu*(1-x**2) ]
45 ts = linspace 1000 (0,50)
46 hi = (ts@>1 - ts@>0)/100
47 sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
48 mplot sol
49
50
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index 2251acd..c243f4b 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -29,7 +29,7 @@ main = mplot (ts : toColumns sol)@
29----------------------------------------------------------------------------- 29-----------------------------------------------------------------------------
30 30
31module Numeric.GSL.ODE ( 31module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..) 32 odeSolve, odeSolveV, ODEMethod(..), Jacobian
33) where 33) where
34 34
35import Data.Packed.Internal 35import Data.Packed.Internal
@@ -41,18 +41,21 @@ import System.IO.Unsafe(unsafePerformIO)
41 41
42------------------------------------------------------------------------- 42-------------------------------------------------------------------------
43 43
44type Jacobian = Double -> Vector Double -> Matrix Double
45
44-- | Stepping functions 46-- | Stepping functions
45data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. 47data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
46 | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'. 48 | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods.
47 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. 49 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
48 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method. 50 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
49 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method. 51 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
50 | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. 52 | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
51 | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points. 53 | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
52 | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. 54 | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
53 | Gear1 -- ^ M=1 implicit Gear method. 55 | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
54 | Gear2 -- ^ M=2 implicit Gear method. 56 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
55 deriving (Enum,Eq,Show,Bounded) 57 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
58
56 59
57-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. 60-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
58odeSolve 61odeSolve
@@ -60,7 +63,7 @@ odeSolve
60 -> [Double] -- ^ initial conditions 63 -> [Double] -- ^ initial conditions
61 -> Vector Double -- ^ desired solution times 64 -> Vector Double -- ^ desired solution times
62 -> Matrix Double -- ^ solution 65 -> Matrix Double -- ^ solution
63odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts 66odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
64 where hi = (ts@>1 - ts@>0)/100 67 where hi = (ts@>1 - ts@>0)/100
65 epsAbs = 1.49012e-08 68 epsAbs = 1.49012e-08
66 epsRel = 1.49012e-08 69 epsRel = 1.49012e-08
@@ -73,11 +76,33 @@ odeSolveV
73 -> Double -- ^ absolute tolerance for the state vector 76 -> Double -- ^ absolute tolerance for the state vector
74 -> Double -- ^ relative tolerance for the state vector 77 -> Double -- ^ relative tolerance for the state vector
75 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) 78 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
79 -> Vector Double -- ^ initial conditions
80 -> Vector Double -- ^ desired solution times
81 -> Matrix Double -- ^ solution
82odeSolveV RK2 = odeSolveV' 0 Nothing
83odeSolveV RK4 = odeSolveV' 1 Nothing
84odeSolveV RKf45 = odeSolveV' 2 Nothing
85odeSolveV RKck = odeSolveV' 3 Nothing
86odeSolveV RK8pd = odeSolveV' 4 Nothing
87odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
88odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
89odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
90odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
91odeSolveV MSAdams = odeSolveV' 9 Nothing
92odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
93
94
95odeSolveV'
96 :: CInt
76 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian 97 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
98 -> Double -- ^ initial step size
99 -> Double -- ^ absolute tolerance for the state vector
100 -> Double -- ^ relative tolerance for the state vector
101 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
77 -> Vector Double -- ^ initial conditions 102 -> Vector Double -- ^ initial conditions
78 -> Vector Double -- ^ desired solution times 103 -> Vector Double -- ^ desired solution times
79 -> Matrix Double -- ^ solution 104 -> Matrix Double -- ^ solution
80odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do 105odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
81 let n = dim xiv 106 let n = dim xiv
82 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) 107 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
83 jp <- case mbjac of 108 jp <- case mbjac of
@@ -86,7 +111,7 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
86 sol <- vec xiv $ \xiv' -> 111 sol <- vec xiv $ \xiv' ->
87 vec (checkTimes ts) $ \ts' -> 112 vec (checkTimes ts) $ \ts' ->
88 createMIO (dim ts) n 113 createMIO (dim ts) n
89 (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' ) 114 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
90 "ode" 115 "ode"
91 freeHaskellFunPtr fp 116 freeHaskellFunPtr fp
92 return sol 117 return sol
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 5c17836..7f1cf68 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -32,7 +32,7 @@
32#include <gsl/gsl_complex_math.h> 32#include <gsl/gsl_complex_math.h>
33#include <gsl/gsl_rng.h> 33#include <gsl/gsl_rng.h>
34#include <gsl/gsl_randist.h> 34#include <gsl/gsl_randist.h>
35#include <gsl/gsl_odeiv.h> 35#include <gsl/gsl_odeiv2.h>
36#include <gsl/gsl_multifit_nlin.h> 36#include <gsl/gsl_multifit_nlin.h>
37#include <string.h> 37#include <string.h>
38#include <stdio.h> 38#include <stdio.h>
@@ -1356,38 +1356,38 @@ int ode(int method, double h, double eps_abs, double eps_rel,
1356 int jac(double, int, const double*, int, int, double*), 1356 int jac(double, int, const double*, int, int, double*),
1357 KRVEC(xi), KRVEC(ts), RMAT(sol)) { 1357 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
1358 1358
1359 const gsl_odeiv_step_type * T; 1359 const gsl_odeiv2_step_type * T;
1360 1360
1361 switch(method) { 1361 switch(method) {
1362 case 0 : {T = gsl_odeiv_step_rk2; break; } 1362 case 0 : {T = gsl_odeiv2_step_rk2; break; }
1363 case 1 : {T = gsl_odeiv_step_rk4; break; } 1363 case 1 : {T = gsl_odeiv2_step_rk4; break; }
1364 case 2 : {T = gsl_odeiv_step_rkf45; break; } 1364 case 2 : {T = gsl_odeiv2_step_rkf45; break; }
1365 case 3 : {T = gsl_odeiv_step_rkck; break; } 1365 case 3 : {T = gsl_odeiv2_step_rkck; break; }
1366 case 4 : {T = gsl_odeiv_step_rk8pd; break; } 1366 case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
1367 case 5 : {T = gsl_odeiv_step_rk2imp; break; } 1367 case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
1368 case 6 : {T = gsl_odeiv_step_rk4imp; break; } 1368 case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
1369 case 7 : {T = gsl_odeiv_step_bsimp; break; } 1369 case 7 : {T = gsl_odeiv2_step_bsimp; break; }
1370 case 8 : {T = gsl_odeiv_step_gear1; break; } 1370 case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
1371 case 9 : {T = gsl_odeiv_step_gear2; break; } 1371 case 9 : {T = gsl_odeiv2_step_msadams; break; }
1372 case 10: {T = gsl_odeiv2_step_msbdf; break; }
1372 default: ERROR(BAD_CODE); 1373 default: ERROR(BAD_CODE);
1373 } 1374 }
1374 1375
1375
1376 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
1377 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
1378 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
1379
1380 Tode P; 1376 Tode P;
1381 P.f = f; 1377 P.f = f;
1382 P.j = jac; 1378 P.j = jac;
1383 P.n = xin; 1379 P.n = xin;
1384 1380
1385 gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; 1381 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
1382
1383 gsl_odeiv2_driver * d =
1384 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
1386 1385
1387 double t = tsp[0]; 1386 double t = tsp[0];
1388 1387
1389 double* y = (double*)calloc(xin,sizeof(double)); 1388 double* y = (double*)calloc(xin,sizeof(double));
1390 int i,j; 1389 int i,j;
1390 int status;
1391 for(i=0; i< xin; i++) { 1391 for(i=0; i< xin; i++) {
1392 y[i] = xip[i]; 1392 y[i] = xip[i];
1393 solp[i] = xip[i]; 1393 solp[i] = xip[i];
@@ -1396,22 +1396,24 @@ int ode(int method, double h, double eps_abs, double eps_rel,
1396 for (i = 1; i < tsn ; i++) 1396 for (i = 1; i < tsn ; i++)
1397 { 1397 {
1398 double ti = tsp[i]; 1398 double ti = tsp[i];
1399 while (t < ti) 1399
1400 { 1400 status = gsl_odeiv2_driver_apply (d, &t, ti, y);
1401 gsl_odeiv_evolve_apply (e, c, s, 1401
1402 &sys, 1402 if (status != GSL_SUCCESS) {
1403 &t, ti, &h, 1403 printf ("error in ode, return value=%d\n", status);
1404 y); 1404 break;
1405 // if (h < hmin) h = hmin; 1405 }
1406 } 1406
1407// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
1408
1407 for(j=0; j<xin; j++) { 1409 for(j=0; j<xin; j++) {
1408 solp[i*xin + j] = y[j]; 1410 solp[i*xin + j] = y[j];
1409 } 1411 }
1410 } 1412 }
1411 1413
1412 free(y); 1414 free(y);
1413 gsl_odeiv_evolve_free (e); 1415 gsl_odeiv2_driver_free (d);
1414 gsl_odeiv_control_free (c); 1416
1415 gsl_odeiv_step_free (s); 1417 return status;
1416 return 0;
1417} 1418}
1419