solve for mathieu duffing type ODE

Kollányi Tibor kollanyi@szgtirix.szgt.uni-miskolc.hu
Thu Mar 20 18:26:00 GMT 2003


Hello,

I would like solve the next mathieu duffing type ODE:
y''+(lambda-2*mu*cos(2*t))y+y^3=0
But there are some problem. Can everybody help me? Thanks.
The code is :

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

int
func (double t, const double y[], double f[],
      void *params)
{
  double mu = *(double *)params;
  double lambda= *((double *)params+1);   
  double epsilon= *((double *)params+2);   
  f[0] = y[1];
  f[1] = -1.0*(lambda-2.0*mu*cos(2.0*t))*y[0]-epsilon*y[0]*y[0]*y[0];
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  double lambda= *((double *)params+1);   
  double epsilon= *((double *)params+2);
  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*(lambda-2.0*mu*cos(2.0*t))-3.0*epsilon*y[0]*y[0]);
  gsl_matrix_set (m, 1, 1, 0);
  dfdt[0] = 0.0;
  dfdt[1] = -4.0*mu*y[0]*sin(2.0*t);
  return GSL_SUCCESS;
}

int
main (void)
{
  const gsl_odeiv_step_type * T
    = gsl_odeiv_step_rk8pd;

  gsl_odeiv_step * s
    = gsl_odeiv_step_alloc (T, 2);
  gsl_odeiv_control * c
    = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e
    = gsl_odeiv_evolve_alloc (2);

  double mu = 0.2;
  double lambda=1.0;
  double epsilon=0.5;
  double *par;
  par=(double*)malloc(3*sizeof(double));
  par[0]=mu;
  par[1]=lambda;   
  par[2]=epsilon;
  gsl_odeiv_system sys = {func, jac, 2, &par};

  double t = 0.0, t1 = 1.0;
  double h = 1e-6;
  double y[2] = { 0.3, 0.1 };

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

      if (status != GSL_SUCCESS)
          break;

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

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);
  return 0;
}
EXECUTING:
.....
----------------------------------------------
1.00000e-06 3.00000e-01 1.00000e-01
6.00000e-06 3.00001e-01 1.00000e-01
3.10000e-05 3.00003e-01 1.00000e-01
1.56000e-04 3.00016e-01 1.00000e-01
7.81000e-04 3.00078e-01 1.00000e-01
7.81000e-04 3.00078e-01 -4.82170e+09
7.81000e-04 3.00078e-01 -8.26019e+09
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan



More information about the Gsl-discuss mailing list