Simple problem with a Runge-Kutta program

This is a discussion on Simple problem with a Runge-Kutta program within the C Programming forums, part of the General Programming Boards category; Hey guys, I am a new member here and I am having a problem with this Runge-Kutta program. I am ...

  1. #1
    Registered User
    Join Date
    Apr 2007
    Posts
    6

    Question 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);
    }

  2. #2
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by civil25 View Post
    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?
    It looks like you've accidentally flip-flopped your RK intermediates in a few of the equations. For instance:

    Code:
    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));
    This should be:

    Code:
    k1[2][j] = delta_t * f1(y1[j]+(k1[1][j]/2));
    k2[2][j] = delta_t * f2(y2[j]+(k2[1][j]/2), y2[j]+(k2[1][j]/2), t[j]+(delta_t/2));
    The same mistake occurs on other lines.
    Last edited by brewbuck; 04-25-2007 at 11:09 AM.

  3. #3
    Registered User
    Join Date
    Apr 2007
    Posts
    6

    Looks right to me

    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.
    Attached Images Attached Images  

  4. #4
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,548
    > print=fopen("Runge-Kutta.xls", "w");
    Native Excel files are not text files, so call it something like .txt

    You could replace all your \t with comma, then output to "Runge-Kutta.csv"
    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.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

  5. #5
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by civil25 View Post
    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.

  6. #6
    Registered User
    Join Date
    Apr 2007
    Posts
    6
    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.

  7. #7
    Registered User
    Join Date
    Apr 2007
    Posts
    6
    What really bothers me is I have a feeling it is going to be soemthing really simple that I keep looking over.

  8. #8
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by civil25 View Post
    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.

  9. #9
    Registered User
    Join Date
    Apr 2007
    Posts
    6
    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.

  10. #10
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by civil25 View Post
    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.
    Holy crap I am a dumb stupid person. The problem is that you multiply by (1/6). This is an integer expression which evaluates to 0!

    You need to multiply by (1.0/6.0)

    EDIT: Fixing that causes the values to change, but they oscillate pretty wildly. You've got some instability for some reason. Also, I call myself dumb because this is the kind of thing I usually see right away. It wasn't until I had it up in the debugger and stepped past the line that it became obvious. Don't feel bad on this one :-)
    Last edited by brewbuck; 04-25-2007 at 11:59 AM.

  11. #11
    Registered User
    Join Date
    Apr 2007
    Posts
    6
    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!

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. problem with a simple program
    By alfredd in forum C++ Programming
    Replies: 4
    Last Post: 02-28-2009, 02:48 PM
  2. Client-server system with input from separate program
    By robot-ic in forum Networking/Device Communication
    Replies: 3
    Last Post: 01-16-2009, 02:30 PM
  3. Running Program Problem
    By warfang in forum C++ Programming
    Replies: 10
    Last Post: 03-28-2007, 02:02 PM
  4. simple login program problem
    By suckss in forum C Programming
    Replies: 11
    Last Post: 11-11-2006, 04:02 PM
  5. Problem with a simple program
    By Salgat in forum C Programming
    Replies: 10
    Last Post: 06-15-2006, 05:57 PM

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