[Top][All Lists]
[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
>