[Help-gsl] Weird behaviour of oderk4

[Alessandro Toso]

Hi all,

I'm writing a solver for systems of ordinary differential equations of
the 2nd order. I use gsl-1.6 and I took the example reported in the
documentation and modified it. I'm interested in simple physical
problems such as mass-spring-damper systems so I wrote my code for a
single equation and it works fine (using rk4 and compared with matlab
routines). The problem comes with systems of coupled equations. The
solutions seem to be good but there is a scale factor equal to the ratio
between the masses of the system.

I attach here my code and I hope that someone could help me.
The code corresponds to these equations:

m_1 x_1" + (k_1+k_2) x_1 - k_2 x_2 = 0
m_2 x_2" - k_2 x_1 + (k_2+k_3) x_2 - k_3 x_3 = 0
m_3 x_3" - k_3 x_2 + k_3 x_3 = 0

Thanks in advance.
Ale
----------------------- Makefile --------------------------
CC = gcc
TARGETS =  oderk4_mimo
         
INCLUDES = /usr/local/include/gsl
CFLAGS = -Wall 

all: $(TARGETS)
        
oderk4_mimo: oderk4_mimo.c
        $(CC) $(CFLAGS) -I$(INCLUDES) -o $@ $< -lm -lgsl -lgslcblas
        
clean:
        rm -fr $(TARGETS)

--------------------- oderk4_mimo.c -----------------------

#include <stdio.h>
#include <sys/time.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

double m[3][3];
double k[3][3];

static int dim = 6;

int func (double t, const double y[], double f[], void *params)
{

  f[0] = - k[0][0] / m[0][0] * y[3] - k[0][1] / m[0][0] * y[4];
  f[1] = - k[1][0] / m[1][1] * y[3] - k[1][1] / m[1][1] * y[4] 
         - k[1][2] / m[1][1] * y[5];
  f[2] = - k[2][1] / m[2][2] * y[4] - k[2][2] / m[2][2] * y[5];

  f[3] = y[0];
  f[4] = y[1];
  f[5] = y[2];

  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, dim, dim);
  
  gsl_matrix * A = &dfdy_mat.matrix; 
    
  gsl_matrix_set (A, 0, 0, 0.0);
  gsl_matrix_set (A, 0, 1, 0.0);
  gsl_matrix_set (A, 0, 2, 1.0);
  gsl_matrix_set (A, 0, 3, 0.0);
  
  gsl_matrix_set (A, 1, 0, 0.0);
  gsl_matrix_set (A, 1, 1, 0.0);
  gsl_matrix_set (A, 1, 2, 0.0);
  gsl_matrix_set (A, 1, 3, 1.0);
  
  gsl_matrix_set (A, 2, 0, -k[0][0] / m[0][0]);
  gsl_matrix_set (A, 2, 1, +k[0][1] / m[0][0]);
  gsl_matrix_set (A, 2, 2, -c[0][0] / m[0][0]);
  gsl_matrix_set (A, 2, 3, -c[0][1] / m[0][0]);
  
  gsl_matrix_set (A, 3, 0, +k[1][0] / m[1][1] * y[0]);
  gsl_matrix_set (A, 3, 1, -k[1][1] / m[1][1]);
  gsl_matrix_set (A, 3, 2, -c[1][0] / m[1][1]);
  gsl_matrix_set (A, 3, 3, -c[1][1] / m[1][1]);
  
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  dfdt[3] = 0.0;
*/
  return GSL_SUCCESS;
}

int main (void)
{

  FILE *fp;
  
  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rk4;

  gsl_odeiv_step * s 
    = gsl_odeiv_step_alloc (T, dim);

  int i = 0;

  double mass[3];
  double stif[3];
  
  fp = fopen("system.dat","r");
  
  fscanf(fp, "%lf", &mass[0]);
  fscanf(fp, "%lf", &mass[1]);
  fscanf(fp, "%lf", &mass[2]);
  fscanf(fp, "%lf", &stif[0]);
  fscanf(fp, "%lf", &stif[1]);      
  fscanf(fp, "%lf", &stif[2]);      
  
  fclose(fp);
  
  m[0][0] = mass[0];
  m[0][1] = 0.0;
  m[0][2] = 0.0;
  
  m[1][0] = 0.0;
  m[1][1] = mass[1];
  m[1][2] = 0.0;
  
  m[2][0] = 0.0;
  m[2][1] = 0.0;
  m[2][2] = mass[2];
 
  k[0][0] =  stif[0] + stif[1];
  k[0][1] = -stif[1];
  k[0][2] =  0.0;
  
  k[1][0] = -stif[1];
  k[1][1] =  stif[1] + stif[2];
  k[1][2] = -stif[2];

  k[2][0] =  0.0;
  k[2][1] = -stif[2];
  k[2][2] =  stif[2];

  gsl_odeiv_system sys = {func, NULL, 4, NULL};

  double t  = 0.0;
  double t1 = 10.0;
  double h  = 1e-3;
  
  double y[dim];
  
  y[0] = 0.0; /* dx1/dt */
  y[1] = 0.0; /* dx2/dt */      
  y[2] = 0.0; /* dx3/dt */      
  y[3] = 1.0; /*  x1    */
  y[4] = 0.0; /*  x2    */
  y[5] = 0.0; /*  x3    */
  
  double y_err[dim];

  double dydt_in[dim];
  
  /* initialise dydt_in from system parameters */
  GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);

  struct timeval t_start, t_end;

  gettimeofday(&t_start, NULL);
  int status = 0;
  
  while (t < t1)
  {
      status = gsl_odeiv_step_apply (s, t, h, 
                                         y, y_err, 
                                         NULL, 
                                         NULL, 
                                         &sys);

      if (status != GSL_SUCCESS)
          break;

      t += h;

      printf ("%f %f %f %f %f %f %f\n", t, y[0], y[1], y[2], y[3], y[4],
y[5]);
      i++;
    }

  gettimeofday(&t_end, NULL);
  
  long time = (t_end.tv_sec * 1e6 + t_end.tv_usec) - 
              (t_start.tv_sec * 1e6 + t_start.tv_usec);
          
  fprintf(stderr,"\n\ntime = %ld us for %d cycles\n\n", time, i);
        
  gsl_odeiv_step_free (s);
  return 0;
}


--------------------- oderk4_mimo.c -----------------------

-- 
Alessandro Toso
Dipartimento di Ingegneria Aerospaziale
Politecnico di Milano
Via La Masa, 34
20156 Milano (Italia)
tel  +39 02 2399 8044
mail address@hidden