help-gsl
[Top][All Lists]

## [Help-gsl] ode problem

 From: Tomasz Samotyjak, PWSZ Subject: [Help-gsl] ode problem Date: Wed, 7 Sep 2005 18:40:52 +0200

```Hi,

I've tried to use gsl fo ode's. The jac and func is:

int
func (double t, const double y[], double f[], void *params)
{
double* mu = (double *)params;
double R1 = mu[0], R2 = R1;
double L1 = mu[1], L2 = L1;
double C1 = mu[2];

f[0] = 1/L1 * ( 400*sin(314*t) - f[1] - R1 * f[0] );
f[1] = 1/C1 * ( f[0] - f[2] );
f[2] = 1/L2 * ( f[1] - R2 * f[2] );

return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
double* mu = (double *)params;
double R1 = mu[0], R2 = R1;
double L1 = mu[1], L2 = L1;
double C1 = mu[2];
gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 3,3);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, -R1/L1);
gsl_matrix_set (m, 0, 1, -1/L1);
gsl_matrix_set (m, 0, 2, 0);

gsl_matrix_set (m, 1, 0, 1/C1);
gsl_matrix_set (m, 1, 1, 0);
gsl_matrix_set (m, 1, 2, -1/C1);

gsl_matrix_set (m, 2, 0, 0);
gsl_matrix_set (m, 2, 1, 1/L2);
gsl_matrix_set (m, 2, 2, -R2/L2);

for( int i=0 ; i<3 ; i++)  dfdt[i] = 0.0;
return GSL_SUCCESS;
}

The main function is:

int
main (void)
{
const gsl_odeiv_step_type * T = gsl_odeiv_step_bsimp;
gsl_odeiv_step * s  = gsl_odeiv_step_alloc (T, 3);
gsl_odeiv_control * c  = gsl_odeiv_control_yp_new (1e-3,0);
gsl_odeiv_evolve * e  = gsl_odeiv_evolve_alloc (3);

double mu[] = {1, 1e-3, 1e-3};
gsl_odeiv_system sys = {func, jac, 3, mu};

double t = 0.0, t1 = 0.20;
double h = 1e-5;
double y[] = { 0.0 , 0.0, 0.0};

ofstream out("wynik.txt");

while (t < t1)
{
int status = gsl_odeiv_evolve_apply (e, c, s,
&sys,
&t, t1,
&h, y);

if (status != GSL_SUCCESS)
break;

cout << t << "\t" << y[0] << "\t" << y[1] << "\t" << y[2] << endl;
out <<  t << "\t" << y[0] << "\t" << y[1] << "\t" << y[2] << endl;
}

gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
out.close();
return 0;
}

But program don't work and prints:

8.192e-015 0 0 0
4.9152e-014 0 0 0
2.53952e-013 0 0 0
1.27795e-012 0 0 0
6.39795e-012 0 0 0
3.1998e-011 0 0 0
1.59998e-010 0 0 0
7.99998e-010 0 0 0
4e-009 -1.#IND -1.#IND -1.#IND
2e-008 -1.#IND -1.#IND -1.#IND
1e-007 -1.#IND -1.#IND -1.#IND
5e-007 -1.#IND -1.#IND -1.#IND
2.5e-006 -1.#IND -1.#IND -1.#IND
1.25e-005 -1.#IND -1.#IND -1.#IND
6.25e-005 -1.#IND -1.#IND -1.#IND
0.0003125 -1.#IND -1.#IND -1.#IND
0.0015625 -1.#IND -1.#IND -1.#IND
0.0078125 -1.#IND -1.#IND -1.#IND
0.0390625 -1.#IND -1.#IND -1.#IND
0.195312 -1.#IND -1.#IND -1.#IND
0.2 -1.#IND -1.#IND -1.#IND

I've tried change abs error and other parameters, but the result is the
same. What to do?

Best
Tomasz Samotyjak

```

reply via email to