Thread: Simple problem with a Runge-Kutta program

Hey guys,

I am a new member here and I am having a problem with this Runge-Kutta program. I am programming in C with the Dev C++ v 4.9.9.2 compiler. For some reason my y1 and y2 values are not changing from 2. Does anyone know what is going on?

insert

Code:

#include <stdio.h>
#include <stdlib.h>

double f1(double y2);

double f2(double y1, double y2, double t);

FILE *print;

int main()
{
 	int j;
	double t_max = 80, delta_t = .2;
	//Define N (size)
	int N = t_max / delta_t;
	double y1[N], y2[N], t[N], k1[5][N], k2[5][N], e[N];
	
	t[0] = 0;
	y1[0] = 2;
	y2[0] = 2;
	
	print=fopen("Runge-Kutta.xls", "w");
	fprintf(print, "Time (s)\ty1\ty2\tk1\tk2\tk3\tk4\terror\n");	
		
	for(j=0; j <= N; j++)
	{        
		k1[1][j] = delta_t * f1(y2[j]);
		k2[1][j] = delta_t * f2(y1[j], y2[j], t[j]);
        			
		k1[2][j] = delta_t * f1(y2[j]+(k2[1][j]/2));
		k2[2][j] = delta_t * f2(y1[j]+(k1[1][j]/2), y2[j]+(k2[1][j]/2), t[j]+(delta_t/2));
		
		k1[3][j] = delta_t * f1(y2[j]+(k2[2][j]/2));
		k2[3][j] = delta_t * f2(y1[j]+(k1[2][j]/2), y2[j]+(k2[2][j]/2), t[j]+(delta_t/2));
		
		k1[4][j] = delta_t * f1(y2[j]+k2[3][j]);
		k2[4][j] = delta_t * f2(y1[j]+k1[3][j], y2[j]+k2[3][j], t[j]+delta_t);
		
		y1[j+1] = y1[j] + (1/6)*(k1[1][j] + 2 * k1[2][j] + 2 * k1[3][j] + k1[4][j]);
		y2[j+1] = y2[j] + (1/6)*(k2[1][j] + 2 * k2[2][j] + 2 * k2[3][j] + k2[4][j]);
		
		t[j+1] = t[j] + delta_t;
        
		fprintf(print, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", t[j], y1[j], y2[j], k2[1][j], k2[2][j], k2[3][j], k2[4][j], e[j]);
	}
	
	fprintf(print, "\n\n");
	fclose(print);
	
    //This line keeps the command prompt open so you have time to read the answer.
    system("PAUSE");

    return 0;
}

//------------------------------------------------------------------------------

double f1(double y2)
{	
	return y2;
}

double f2(double y1, double y2, double t)
{
	double c = .75, m = 40, k = .7, f = 2, w = .2;
	
    return (-c/m)*y2-(k/m)*y1+((f*sin(w*t))/m);
}

I just checked it against a working Mathcad program and it looks right to me, here is an image of what is in Mathcad because I might have missed something. Thank you for the help, but I still don't understand what is wrong with it.

Originally Posted by civil25

I just checked it against a working Mathcad program and it looks right to me, here is an image of what is in Mathcad because I might have missed something. Thank you for the help, but I still don't understand what is wrong with it.

Dangit, I can't see what's wrong. You're right, I misread that mess of expressions. If this is NOT for an assignment I can show you some working code.

This was for an assignment, I already turned it in, it was due last Monday. Now this is just out of aggrevation b/c I don't like doing something and getting this close but not finishing it.

What really bothers me is I have a feeling it is going to be soemthing really simple that I keep looking over.

Originally Posted by civil25

What really bothers me is I have a feeling it is going to be soemthing really simple that I keep looking over.

One thing I notice is that your expressions for computing the k's are all multiplied by delta_t on the outside. I've never seen it done that way. For instance look at the Wikipedia page on this topic: http://en.wikipedia.org/wiki/Runge-Kutta_methods

There is no outside multiplication by delta_t in the k formulas. I've never seen that before.

It is the same as what I did only I multiplied delta_t in each of the k formulas instead of doing it only once in the y formula. Granted the way that is presented on wikipedia is less time consuming, they will give you the same exact answer.

That did it! But now I just have to figure out why it is off from the analitical solution... But that will have to wait until after my finals... Thank you very much, atleast now I can work with it!