# 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

```