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

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

Join Date
Dec 2012
Posts
1
Rep Power
0

#### 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:

Code:
```#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);

main()

{
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");
fclose(output);
}

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. Please use [code][/code] tags when posting code.
If you're quick, you might still be able to edit your first post.
3. No Profile Picture
Still Learning
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Location