[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] odeiv2 rk2imp driver time step

From: Tuomo Keskitalo
Subject: Re: [Help-gsl] odeiv2 rk2imp driver time step
Date: Sun, 01 Jan 2012 17:34:28 +0200
User-agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv: Gecko/20111108 Thunderbird/3.1.16


thanks for the example code! It seems that the modified Newton iteration method of rk*imp (modnewton1.c) starts to fail at t=14.875 (iteration does not improve the solution). Since both bsimp and msbdf seem to work OK, I think that modnewton1 just is not appropriate for your problem. We don't have another non-linear equation solver implemented for rk*imp steppers, so I suggest you use either bsimp or msbdf, if you can. Generally speaking, you can compare results from different steppers and error tolerances to get some confidence that your results are valid.

odeiv2 functions can make multiple user function calls at a time point. E.g. non-embedded Runge-Kutta steppers use step size halving for error estimation on each step. Those function calls are necessary.


On 12/22/2011 09:55 PM, Farkas, Illes wrote:

Hi Juan,

Sorry for the delay and thanks for all the help. I've just finished
trimming it down to a simple test case. Please see the attached source
(test.c), the stdout (test.out.txt), the gnuplot cmd file (test.gnu) and
the resulting plot (, test.pdf). It is very likely that I'm making
some very dumb mistake. I just can't find it. Thanks for any help !

The compilation command was: gcc test.c -Wall -O3 -I...<include>...
-L...<lib>... -lm -lgsl -lgslcblas -o test

The test output file (test.out.txt) contains output from the main program
(lines starting with "main") and output from the function (lines starting
with "ode_system_rhs_function").

According to the last ~50 lines of test.out.txt it seems that the attempted
time step is halved many times.

One more question. According to test.out.txt the r.h.s. function is
accessed 3 or 2 times at each time point. Can this be reduced to accessing
only once (in order to speed up rk2imp) ?


2011/12/22 Juan Pablo Amorocho D.<address@hidden>


I'm curious how you solved your problem. I would be grateful if you post
your solution on the help-gsl list.

-- Juan

2011/12/21 Farkas, Illes<address@hidden>

Thanks, Tuomo

I just ran into something unexpected. I am trying to find out where
I'm making a mistake.

I'm integrating with rk2imp a 3d ODE that has constants, linear and second
order polynomial terms on the r.h.s. I use gsl_odeiv2_driver_apply to
evolve the ODE in steps of 0.125s. After 10s or so (the exact time varies
with the parameters) all three variables converge (very little relative
change). However, some time (>  5s ) later there is an update
when gsl_odeiv2_driver_apply returns the FAILURE value: -1. After logging
the current time directly from the "function" and "jacobian" (used by
the gsl_odeiv2_system, which is driver by the driver), I found that this
particular update fails, because the time step is halved again and again
until it reaches the limit of numerical precision.

Have you seen a similar error before ?

I have the Ascher Petzold book. Will sections 5.4.3 (modified Newton
iteration) and 4.7 (implicit methods) be helpful for this problem or shall
I use a different resource ?


2011/12/1 Tuomo Keskitalo<address@hidden>


rk2imp (among other implicit methods) in ode-initval2 uses a modified
Newton iteration instead of old functional iteration in solving the
of non-linear equations. E.g. for stiff problems Newton iteration is
more powerful. Because of that, the ODE-solver can use larger step
Newton iteration converges with larger step sizes, while functional
iteration does not.

On 12/01/2011 08:42 PM, Farkas, Illes wrote:


I just tested the speed of rk2imp with the simple harmonic oscillator
(dx/dt=-y, dy/dt=x). In the first test I used
(does *not* use Jacobian) and in the 2nd test I used
simply gsl_odeiv2_driver (uses Jacobian). With the same parameters the
first version ran for 26s and the second version finished below 1s. Is
test wrong? Or is there really such a big difference?

Help-gsl mailing list


Help-gsl mailing list

Help-gsl mailing list


reply via email to

[Prev in Thread] Current Thread [Next in Thread]