Forums: » Register « |  Free Tools |  User CP |  Games |  Calendar |  Members |  FAQs |  Sitemap |  Support |

New Free Tools on Dev Shed!

#1
December 7th, 2012, 01:08 PM
 owlman76
Registered User

Join Date: Dec 2012
Posts: 1
Time spent in forums: 27 m 38 sec
Reputation 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
December 7th, 2012, 02:25 PM
 salem
Contributed User

Join Date: Jun 2005
Posts: 4,252
Time spent in forums: 2 Months 4 Weeks 1 Day 12 h 25 m 37 sec
Reputation Power: 1809
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

#3
December 16th, 2012, 06:00 AM
Still Learning

Join Date: Dec 2012
Posts: 55
Time spent in forums: 18 h 57 m 51 sec
Reputation Power: 38
It's been done.

This has probably been solve before. A goggle search turned up
http://people.sc.fsu.edu/~jburkardt/c_src/rk4/rk4.html

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

 Viewing: Dev Shed Forums > Programming Languages > C Programming > Need help modifying my C program to calculate error in rk4 method!!