help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] ODE example, passing parameters


From: gottstew
Subject: [Help-gsl] ODE example, passing parameters
Date: Thu, 10 Jun 2010 17:02:00 +0200
User-agent: SquirrelMail/1.4.19

Dear Sir or Madam,

I slightly modified your ODE example and tried to pass multiple
parameters. As you recommended I defined a structure in order to do that.
But now I always get a warning in the following line:
gsl_odeiv_system sys = {func, jac, 3, pz};

where pz is a pointer to the structure that contains multiple paramters:

"warning: initialization from incompatible pointer type"

I get this warning twice. I guess the first time when I call the function
"func" and the second time when I call "jac".

Here is the code I use (is also attached):

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




struct par{
        double mu;
        double om;
        double de;
};

int
func (double t, const double y[], double f[], struct par *pz)
{

  double mu=pz->mu;
  double om=pz->om;
  double de=pz->de;
  //printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);

  f[0] = om*y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  f[2] = de*y[1]*y[2];

  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
double dfdt[], struct par *pz)
{
  double mu=pz->mu;
  double om=pz->om;
  double de=pz->de;
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 3, 3);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, om);
  gsl_matrix_set (m, 0, 2, 0.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  gsl_matrix_set (m, 1, 2, 0.0);
  gsl_matrix_set (m, 2, 0, 0.0);
  gsl_matrix_set (m, 2, 1, de*y[2]);
  gsl_matrix_set (m, 2, 2, de*y[1]);

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  return GSL_SUCCESS;
}
/*
void print_para(struct par *pz)
{
        printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);
}
*/

int
main (void)
{
  FILE *fp;
  fp = fopen("odetest.dat","w");
  if(fp == NULL)
   {
      fprintf(stderr,"error: can't open odetest.dat!\n");
      exit(EXIT_FAILURE);
   }



  const gsl_odeiv_step_type * T
    = gsl_odeiv_step_rk8pd;

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


  struct par z, *pz;
  pz=&z;
  z.mu = 10;
  pz->om = 1.45;
  pz->de = 0.2;

  printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);

  //print_para(pz);

  gsl_odeiv_system sys = {func, jac, 3, pz};

  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y[3] = { 2.0, 2.0, 2.0 }; //initial cond

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

         if (status != GSL_SUCCESS)
                break;

         fprintf (fp,"%.5e %.5e %.5e %.5e\n", t, y[0], y[1], y[2]);

    }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);


  return 0;
}


Thanks in advance!

Best,
Willi





Attachment: test_exa2.c
Description: Text Data


reply via email to

[Prev in Thread] Current Thread [Next in Thread]