GSL ode-initval development

GSL ode-initval development

[Tuomo Keskitalo]

Hello all,

I've been thinking about adding some implicit one step solvers for GSL 
ode-initval. However, I think the current ode-initval framework is not 
too easy to use for implicit solvers.

Implicit solvers result in a group of algebraic non-linear equations, 
for which there exists several solution strategies. The choice of the 
efficient strategy depends on the problem. For example, some modified 
Newton iteration methods are suitable for stiff systems, while 
functional iteration works for non-stiff systems. I would like to give 
the user the freedom to choose the non-linear eq solver separately from 
the stepping method. Does anyone see a way to do this with current 
framework?

I am currently considering to add a new framework part for specifying 
the non-linear equation solver. However, this would break the current 
framework. I wonder if I should write a separate "ode-initval2" library 
because of this. Any comments?

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: GSL ode-initval development

[Brian Gough]

At Sun, 17 Aug 2008 11:38:15 +0300,
Tuomo Keskitalo wrote:
> Implicit solvers result in a group of algebraic non-linear
> equations, for which there exists several solution strategies. The
> choice of the efficient strategy depends on the problem. For
> example, some modified Newton iteration methods are suitable for
> stiff systems, while functional iteration works for non-stiff
> systems. I would like to give the user the freedom to choose the
> non-linear eq solver separately from the stepping method. Does
> anyone see a way to do this with current framework?

If it is only a small number (and they don't need any additional
parameters) then it's simplest to use the existing framework with
explicit names for each combination, e.g. gsl_odeiv_step_foo_newt,
gsl_odeiv_step_foo_imp, etc

We mainly need additional stiff solvers, so in practical terms the
functional iteration case probably isn't needed -- the non-stiff case
should be covered adequately by the existing RK-type rules.

Of course, in terms of internal implementation it's better to have the
methods interchangeable and it is worth implementing it that way if
it's not too complicated -- but in terms of exposing them to the user
I don't think it's necessary to have all the permutations.

> I am currently considering to add a new framework part for specifying
> the non-linear equation solver. However, this would break the
> current framework. I wonder if I should write a separate
> "ode-initval2" library because of this. Any comments?

In terms of testing and being useful in existing code, it's a lot
easier to go with the current framework.

-- 
Brian Gough

Re: GSL ode-initval development

[Tuomo Keskitalo]

On 08/18/2008 09:09 PM, Brian Gough wrote:

> At Sun, 17 Aug 2008 11:38:15 +0300,
> Tuomo Keskitalo wrote:
>> Implicit solvers result in a group of algebraic non-linear
>> equations, for which there exists several solution strategies.
> If it is only a small number (and they don't need any additional
> parameters) then it's simplest to use the existing framework with
> explicit names for each combination, e.g. gsl_odeiv_step_foo_newt,
> gsl_odeiv_step_foo_imp, etc

I think there will not be too many (hopefully), and your suggestion is OK.

I have been writing a non-linear equation solver based on a modified 
Newton iteration method, and a question about failures has come up. What 
value should the non-linear equation solver return to the stepper (and 
from there to user) in the case that the method reaches a maximum number 
of Newton iterations and therefore fails? GSL_CONTINUE or GSL_FAILURE or 
something else?

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

ode-initval stepper logic

[Tuomo Keskitalo]

Hello,

which values of variables y should a stepper which uses step doubling 
for error estimation return: The values from the single step or the 
values from two half steps? You would expect to get the values of single 
step from "a stepper", but the values from two half steps should be more 
accurate, and they are available.

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: GSL ode-initval development

[Brian Gough]

At Sun, 31 Aug 2008 10:25:26 +0300,
Tuomo Keskitalo wrote:
> I have been writing a non-linear equation solver based on a modified 
> Newton iteration method, and a question about failures has come up. What 
> value should the non-linear equation solver return to the stepper (and 
> from there to user) in the case that the method reaches a maximum number 
> of Newton iterations and therefore fails? GSL_CONTINUE or GSL_FAILURE or 
> something else?

I'd suggest GSL_FAILURE to start with. Hopefully the solver won't fail
that often.

> which values of variables y should a stepper which uses step doubling 
> for error estimation return: The values from the single step or the 
> values from two half steps? You would expect to get the values of single 
> step from "a stepper", but the values from two half steps should be more 
> accurate, and they are available.

The final value of y should come from the two half steps as they are
more accurate.  See rk4.c for an example and the error estimate.

-- 
Brian Gough