1. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Dec 2012
    Rep Power

    Need help modifying my C program to calculate error in rk4 method!!

    Hello everyone
    I have constructed my program to perform runge-kutta 4th order method to solve a simple harmonic oscillator. I need to construct the program though to calculate the exact solution as well and calculate the global error.

    Here is my program so far:

    #include <stdio.h>
    #include <math.h>
    #include <stdlib.h>
    #define N 2                                   
    #define MIN 0.0
    #define K 1            /* spring constant */
    #define M 1            /* mass of block */
    void runge4(double x, double y[], double step);
    double f(double x, double y[], int i);
       double x, y[N], v0, x0, tmax, dt;
       int j, n;
       FILE *output;                             /* save data in pa2.dat */
       output = fopen("pa2.dat","w");
       printf("\nPlease Enter initial position and velocity X0 and V0: ");
       scanf("%le %le", &x0, &v0);
       printf("\nPlease Enter Tmax and N: ");
       scanf("%le %d", &tmax, &n);
       dt = tmax/(n-1);
       y[0] = x0;                                /* initial position  */
       y[1] = v0;                                /* initial velocity  */
       fprintf(output, "%f\t%f\n", x, y[0]);
       for(x = MIN; x <= tmax ; x += dt)
          runge4(x, y, dt);
          fprintf(output, "%f\t%f\n", x, y[0]);   /* position vs. time */
       printf("data stored in rk4.dat\n\n");
    void runge4(double x, double y[], double step)
       double h=step/2.0,                        
              t1[N], t2[N], t3[N],               
              k1[N], k2[N], k3[N],k4[N];
       int i;
       for (i=0; i<N; i++) t1[i] = y[i]+0.5*(k1[i]=step*f(x, y, i));
       for (i=0; i<N; i++) t2[i] = y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
       for (i=0; i<N; i++) t3[i] = y[i]+    (k3[i]=step*f(x+h, t2, i));
       for (i=0; i<N; i++) k4[i] =                 step*f(x + step, t3, i);
    /* equations defined */
    double  f(double x, double y[], int i)
       if (i == 0) return((2*y[0]*(1 - y[1]));               /* RHS of first equation */
       if (i == 1) return((-y[1]*(1 - y[0]));        /* RHS of second equation */
    The exact solutions for the function are:

    x(t) = x0cos(wt) + v0/w*sin(wt)
    v(t) = -wx0*sin(wt) + v0*cos(wt)

    I obviously know my program lacks true time dependence which is needed but I am pretty clueless on how to actually include the exact solutions and calculate the global error between the two. Any pointers would be appreciated.
  2. #2
  3. Contributed User
    Devshed Specialist (4000 - 4499 posts)

    Join Date
    Jun 2005
    Rep Power
    Please use [code][/code] tags when posting code.
    If you're quick, you might still be able to edit your first post.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper
  4. #3
  5. No Profile Picture
    Still Learning
    Devshed Newbie (0 - 499 posts)

    Join Date
    Dec 2012
    Montreal, Canada
    Rep Power

    It's been done.

    This has probably been solve before. A goggle search turned up

    You may have to look at a few other samples to get exactly what you are looking.

IMN logo majestic logo threadwatch logo seochat tools logo