help-gsl
[Top][All Lists]
Advanced

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

Re: [Help-gsl] convergence check in odeiv_rk4imp ?


From: Forest Yang
Subject: Re: [Help-gsl] convergence check in odeiv_rk4imp ?
Date: Wed, 6 Jan 2010 08:50:57 -0500

Hi,

   I did also use the driver_apply function, which I guess is not
fixed stepsize from the signature of the function
"odeiv2_driver_apply(d, &t, ti, y);"
step_apply looks more like the fixed stepsize, as shown in odeiv2's example.

   I will look at the code see what happens there.

  Thanks.


Forest



int
func (double t, const double y[], double f[],
      void *params)
{
    printf("#t= %.5e\n", t);
  f[0] = y[1];
  f[1] = -y[0];
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -1.0);
  gsl_matrix_set (m, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}


int main (void)
{
    int i = 0;
    const odeiv2_step_type * T = odeiv2_step_rk4imp;

    odeiv2_step * s = odeiv2_step_alloc (T, 2);

    double mu = 10;
    odeiv2_system sys = {func, jac, 2, &mu};

    double t = 0.0, t1 = 2.0*M_PI;
    double h = M_PI*2.0/3.0;
    double y[2] = { 0.0, 1.0 }, y_err[2];
    double dydt_in[2], dydt_out[2];

    odeiv2_driver *d = odeiv2_driver_alloc_y_new(&sys,
odeiv2_step_rk4imp, 10, 10, 0.0);
    odeiv2_driver_set_hmin(d, h);
    odeiv2_step_set_driver(s, d);

    FILE* fp = fopen("arnold_cat_sup.dat", "r");
    int dtmp, idx;
    while (!feof(fp)) {
        fscanf(fp, "%lf %lf %d", &y[0], &y[1], &dtmp);
        t = 0;
    /* initialise dydt_in from system parameters */
    GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);

    printf ("%.5e %.5e %.5e %d\n", t, y[0], y[1], idx);
    //for (i = 1; i <= 3; ++i) {
    while (t < t1) {
        double ti = i*t1/3;
        int status = odeiv2_step_apply (s, t, h, y, y_err, dydt_in,
dydt_out, &sys);
        //int status = odeiv2_driver_apply(d, &t, ti, y);

        if (status != GSL_SUCCESS)
            break;

        dydt_in[0] = dydt_out[0];
        dydt_in[1] = dydt_out[1];

        t += h;

        printf ("%.5e %.5e %.5e %d\n", t, y[0], y[1], idx);
    }
    ++idx;
    }

    odeiv2_step_free (s);
    return 0;
}


On Wed, Jan 6, 2010 at 7:15 AM, Tuomo Keskitalo <address@hidden> wrote:
> Hello,
>
> My guess is that rk4imp can't make Newton iteration converge when you use a
> large step size. You should be able to use a small fixed step size with any
> stepper, if that is really needed.
>
> If fixed step sizes are not mandatory, I suggest you use the driver or
> evolve functions to traverse from t0 to t1. Please see the example codes in
> the documentation.
>
> Regards,
> Tuomo
>
> On 01/06/2010 01:34 AM, Forest Yang wrote:
>>
>> Hi
>>
>>  I was trying to evaluate the rk4imp solver at large step, but not
>> being able to do it on a fixed stepsize way.
>> Can I truely do rk4imp with a fixed time step ?
>>
>> I used the following two lines, and also put a print line in f, the
>> time printed out have many steps between [t,t+h], Is it normal in the
>> internal stage when solving the fixed point or still an adaptive
>> stepsize control ?
>>
>>
>>  GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);
>>  int status = odeiv2_step_apply (s, t, h, y, y_err, dydt_in, dydt_out,
>> &sys);
>>
>> Lingyun
>>
>>
>>
>> On Thu, Dec 31, 2009 at 1:24 AM, Tuomo Keskitalo <address@hidden>
>> wrote:
>>>
>>> Hello,
>>>
>>> take a look at my ode-initval2 extension to GSL, it's rk4imp iterates to
>>> convergence. Note: I'm currently working on finetuning the ode-initval2
>>> framework, it will soon change slightly from version 0.9.
>>>
>>> http://iki.fi/tuomo.keskitalo/gsl/ode-initval2/
>>>
>>> On 12/29/2009 08:35 PM, Forest Yang wrote:
>>>
>>>> Hi All,
>>>>
>>>>   Is there any plan for implementing convergence check in implicit
>>>> RK4 method ?
>>>> It seems the 2 stages Gauss points are only iterate 3 steps, without
>>>> checking whether it is converged or not.
>>>> Another possible way could be using the Jacobian to solve for "Y_i",
>>>> the fixed point.
>>>>
>>>>  The reason I care about this is, Implicit Gauss RK4 is symplectic,
>>>> which for Hamiltonian Dynamics simulation is a very important
>>>> property,
>>>> it has better conservation of the first integral.
>>>>
>>>>
>>>> Best regards,
>>>> Forest.
>>>>
>>>>
>>>> _______________________________________________
>>>> Help-gsl mailing list
>>>> address@hidden
>>>> http://lists.gnu.org/mailman/listinfo/help-gsl
>>>>
>>>
>>> --
>>> address@hidden
>>> http://iki.fi/tuomo.keskitalo
>>>
>>
>
> --
> address@hidden
> http://iki.fi/tuomo.keskitalo
>




reply via email to

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