Like Tree1Likes
  • 1 Post By tabstop

Numerical integration using rkf45

This is a discussion on Numerical integration using rkf45 within the C Programming forums, part of the General Programming Boards category; Hi all, I am writing a C program to numerically integrate a differential equation. This differential equation corresponds to dynamics ...

  1. #1
    Registered User
    Join Date
    Dec 2010
    Posts
    16

    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 ???

  2. #2
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Why do you think rkf45 requires the Jacobian?

  3. #3
    Registered User
    Join Date
    Dec 2010
    Posts
    16
    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;
    }

  4. #4
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    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.
    Salem likes this.

  5. #5
    Registered User
    Join Date
    Dec 2010
    Posts
    16
    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.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Replies: 2
    Last Post: 10-27-2009, 08:19 AM
  2. Integration with Web?
    By dluthcke in forum C++ Programming
    Replies: 25
    Last Post: 09-30-2009, 12:02 AM
  3. Numerical Integration by Trap Rule.
    By hexagons in forum C Programming
    Replies: 3
    Last Post: 03-25-2009, 10:12 AM
  4. Numerical Integration
    By Crazed in forum C Programming
    Replies: 13
    Last Post: 03-03-2008, 02:01 PM
  5. Integration
    By westclok in forum C Programming
    Replies: 8
    Last Post: 02-01-2005, 04:06 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21