Thread: Numerical integration using rkf45

Hi all,
I am writing a C program to numerically integrate a differential equation. This differential equation corresponds to dynamics of a system. I am using rkf45 to integrate it. But the problem with rkf45 is that it requires Jacobian of the equations. I want to get rid of jacobian. ode45 in MATLAB does the same thing w/o using Jacobian. Is there any way to get rid of jacobian in C ???

Why do you think rkf45 requires the Jacobian?

This is the code I am using to integrate a function indicated by "func_aug" and its jacobian is indicated by "jac_aug". Thus jacobian is required as an argument in the following integrator. Can you suggest me how to get rid of it??

Code:

void ode45(double tStart,double tEnd,int dim,double *x0,int *steps)
{
   const gsl_odeiv_step_type *T = gsl_odeiv_step_rkf45;
   gsl_odeiv_step *s = gsl_odeiv_step_alloc(T,3);
   gsl_odeiv_control *c = gsl_odeiv_control_y_new(1e-10,0.0);
   gsl_odeiv_evolve *e  = gsl_odeiv_evolve_alloc(dim);
   gsl_odeiv_system sys = {func_aug,jac_aug,dim,steps};
  
      
   double t = tStart, t1 = tEnd;
   double h = 1e-9;
   int ctr = 0;
   while(t<t1)
   {
      int status = gsl_odeiv_evolve_apply(e,c,s,&sys,&t,t1,&h,x0);
      if(status != GSL_SUCCESS)
      {
         printf("\nFailure from evolve!!!!!!!!\n");
         break;
      }
      ctr++;
   }
   *steps = ctr-1; 
   gsl_odeiv_evolve_free(e);
   gsl_odeiv_control_free(c);
   gsl_odeiv_step_free(s);
   return;
}

So: hopefully you recognize the names as coming from the GNU Scientific Library. Have you looked in the GSL manual? Specifically 26.1 where it says

int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
This function should store the vector of derivative elements in the array dfdt and the Jacobian matrix J_ij in the array dfdy, regarded as a row-ordered matrix J(i,j) = dfdy[i * dimension + j] where dimension is the dimension of the system.

Not all of the stepper algorithms of gsl_odeiv2 make use of the Jacobian matrix, so it may not be necessary to provide this function (the jacobian element of the struct can be replaced by a null pointer for those algorithms).

(emphasis added)
26.2 confirms that rkf45 is one of those algorithms.

Thanks a lot. I was really trying to get rid of this jacobian as the equations are very complex and the order of the system is high. I found that "gsl_odeiv2_step_apply' does not explicitly require jacobian.