Runge-Kutta Problem

nickbrown05 · 11-13-2008, 07:45 AM

Project
Use the fourth order Runge-Kutta algorithm to solve the differential equation.

dy/dx = -y, y(0) = 1

thats the problem baiscally, below is the code I have got so far and so far as I am a complete beginner to c/c++ I'm having great difficulty getting this to work.
Please any help on this would be greatly appreciated as i currently feeling like im banging my head against a large c shape wall

Code:

#include<stdio.h>
#define xs 0.1
#define xf 5
FILE *output; 

main()
{
    double t, y, h, y0 , yn;
    double f(double x, double y);
    double runge(double x, double y)
    int j;
output=fopen("data3.dat", "w"); /* external filename */
    y0 = 1;
    h = 0.1;
    yn = y0;
    fprintf(output, "0\t%f\n", y);

    for (j=0; j*xs<=xf; j++)             /* the time loop */
    {
       t=j*xs;
          y-=runge(t, y,);

     fprintf(output, "%f\t%f\n", t, y);
     }

fclose(output);
}
   
double runge(double x, double y)
{
    double K1    = (H * f(x,y));
    double K2    = (H * f((x + 1 / 2 * H), (y + 1 / 2 * K1)));
    double K3    = (H * f((x + 1 / 2 * H), (y + 1 / 2 * K2)));
    double K4    = (H * f((x + H), (y + K3)));
	double runge = (y + (1 / 6) * (K1 + 2 * K2 + 2 * K3 + K4));
	return runge;
}
   double f(double x, double y)
{
	double f = (-y);
	return f;
}

matsp · 11-13-2008, 07:53 AM

I can see some things that seem wrong, but perhaps you can explain what it is that isn't working, so we can fix THAT, rather than something that you may well already be on the way of fixing (or not important at this point).

My guess is that it doesn't compile because H is unknown (I haven't even tried to compile it).

It's ALLOWED to use the same name for a function and a local variable in the function, but it's a bad idea. Rename your "runge" variable inside the runge function to something else (result, perhaps?)

--
Mats

nickbrown05 · 11-13-2008, 08:24 AM

Ok my goal it to creat a table for all the values of y so that I can then take the data values and them plot a graph, i have attempted the code again beneath should be a better code then my first but still having major troubles...
Half the table I produce of x from0 to 5 at increments of 0.1 that comes up fine.... but I should get alterning values for y but keep just getting -1 and so obviously the functons input arent running.... so again any help with the code would be awesome

when compiled im greeted with: (for the beneath code and havent a clue where i should go from here)

39 [Warning] previous implicit declaration of `f'
48 type mismatch with previous implicit declaration
48 [Warning] `f' was previously implicitly declared to return `int'

Code:

//* A Runge-Kutta Method for solving Differential Equations*/
/* dy/dt = -y(t), y(0)=1, 0<=t<=5, start h=0.1*/

#include <stdio.h>
 
#define dist 0.1		/* stepsize in t */
#define MAX 5			/* max for t */ 
 
FILE *output;			/* internal filename */

main()
{
    double t, y;
    double rkutta( double x, double y, double h);   /*Runge Kutta Function */
    double f(double x, double y);                   /*Function derivative*/
    int n;

    output=fopen("ODE3.dat", "w");	/* External filename */

    y=1;				            /* Initial condition */
    fprintf(output, "0\t%f\n", y);

    for (n=0;dist*n<=MAX;n++)	    /* The time loop */
    {
       t=n*dist;
          y-=rkutta(t, y, dist);
 
   fprintf (output, "%f\t%f\n", t, y);
   }

   fclose(output);
}                   /* End of main function*/

double rkutta(double x, double y, double h)  /*Called on RK function*/
{
    double k1, k2, k3, k4;
    double H = h/2;
    
    k1 = (h*f(x,y));
    k2 = (h*f(x+H,y+(k1/2)));
    k3 = (h*f(x+H,y+(k2/2)));
    k4 = (h*f(x+H,y+k3));
    
    return(y+(k1+2*k2+2*k3+k4)/6);
}

double f(double x, double y)
{
return(-y);
}

tabstop · 11-13-2008, 09:22 AM

Function prototypes need to appear above the declaration of main, not inside main. Also, as I remember RK4, k4 needs to be evaluated at x+h, not x+H. Also, you have written rk to return the new value of y, but you subtract it from y in your loop instead of replacing the value.

nickbrown05 · 11-13-2008, 09:32 AM

ok how would i go about doing that then???????????

tabstop · 11-13-2008, 09:35 AM

Well.... put your function prototypes above main, not inside it; change x+H to x+h in your evaluation of k4; and replace y with the function call value, not subtract the function call value from y.

nickbrown05 · 11-13-2008, 09:53 AM

thank you for your reply, but because im such a noob i dont quite get what you're saying with "replace y with the function call value, not subtract the function call value from y."........just realised may off down off sounding abit like an a-hole previously sorry

tabstop · 11-13-2008, 10:28 AM

The line

Code:

y-=rkutta(t, y, dist);

takes the function call and subtracts it from y. You don't want to subtract it from y, you want it to replace y, so don't subtract it.

nickbrown05 · 11-13-2008, 10:49 AM

This is how im currently looking so thankyou as the program is working better then before but im still having troubles with y i should be producing a results that when plotted look like an exp decrease , but at least now they're plotting something any further suggestions would result in a prize of me hailing you "King or Queen" depending on your mood

Code:

/* A Runge-Kutta Method for solving Differential Equations*/
/* dy/dt = -y(t), y(0)=1, 0<=x<=5, start h=0.1*/

#include <stdio.h>

#define dist 0.1		/* stepsize */
#define xf 5			/* max for x */ 

 
FILE *output;			/* internal filename */

double rkutta( double x, double y, double h);   /*Runge Kutta Function */
double F(double x, double y);                   /*Function derivative*/

main()
{
    double t, y, h;
    int n;
  
    output=fopen("xydata.dat", "w");	/* External filename */
    h=0.1;
    y=1;				            /* Initial condition */
    fprintf(output, "0\t%f\n", y);

    for (n=0;dist*n<=xf;n++)	    /* The time loop */
    {
       t=n*dist;
          y=rkutta(t, y, dist);
 
   fprintf (output, "%f\t%f\n", t, y);
   }

   fclose(output);
}                   /* End of main function*/

double rkutta(double x, double y, double h)  /*Called on RK function*/
{
    double yn, k1, k2, k3, k4;
    double H = h/2.0;
          
    k1 = (h*F(x, y));
    k2 = (h*F(x+H, y+(k1/2)));
    k3 = (h*F(x+H, y+(k2/2)));
    k4 = (h*F(x+h, y+k3));
    return(y+=(y+(k1+2*k2+2*k3+k4)*1/6));
}
   

double F(double x, double y)                /*Called on derivative*/
{
    if (x==0){
    return(-y);
}
}

tabstop · 11-13-2008, 10:53 AM

So now you're adding y twice in the RK function itself. You should only add it once (which is what you had when you started). Also F should return -y no matter what x is -- right now you're returning "?" most of the time from F, and who knows what's going to happen when you start adding "?" to things.

nickbrown05 · 11-13-2008, 11:00 AM

tabstop EST KING KING KING.............................................. .........................................

IT EFFING WORKS, for the first time I can see a bond between man and machine no longer will I live in fear of being snack attacked by the terminator

11-13-2008, 07:53 AM	#2
matsp Kernel hacker Join Date: Jul 2007 Location: Farncombe, Surrey, England Posts: 15,686	I can see some things that seem wrong, but perhaps you can explain what it is that isn't working, so we can fix THAT, rather than something that you may well already be on the way of fixing (or not important at this point). My guess is that it doesn't compile because H is unknown (I haven't even tried to compile it). It's ALLOWED to use the same name for a function and a local variable in the function, but it's a bad idea. Rename your "runge" variable inside the runge function to something else (result, perhaps?) -- Mats __________________ Compilers can produce warnings - make the compiler programmers happy: Use them! Please don't PM me for help - and no, I don't do help over instant messengers.
matsp is offline

11-13-2008, 09:22 AM	#4
tabstop and the Hat of Guessing Join Date: Nov 2007 Posts: 8,740	Function prototypes need to appear above the declaration of main, not inside main. Also, as I remember RK4, k4 needs to be evaluated at x+h, not x+H. Also, you have written rk to return the new value of y, but you subtract it from y in your loop instead of replacing the value.
tabstop is offline

11-13-2008, 09:32 AM	#5
nickbrown05 Registered User Join Date: Nov 2008 Posts: 6	ok how would i go about doing that then???????????
nickbrown05 is offline

11-13-2008, 09:35 AM	#6
tabstop and the Hat of Guessing Join Date: Nov 2007 Posts: 8,740	Well.... put your function prototypes above main, not inside it; change x+H to x+h in your evaluation of k4; and replace y with the function call value, not subtract the function call value from y.
tabstop is offline

11-13-2008, 09:53 AM	#7
nickbrown05 Registered User Join Date: Nov 2008 Posts: 6	thank you for your reply, but because im such a noob i dont quite get what you're saying with "replace y with the function call value, not subtract the function call value from y."........just realised may off down off sounding abit like an a-hole previously sorry
nickbrown05 is offline

11-13-2008, 10:28 AM	#8
tabstop and the Hat of Guessing Join Date: Nov 2007 Posts: 8,740	The line Code: y-=rkutta(t, y, dist); takes the function call and subtracts it from y. You don't want to subtract it from y, you want it to replace y, so don't subtract it.
tabstop is offline

11-13-2008, 10:53 AM	#10
tabstop and the Hat of Guessing Join Date: Nov 2007 Posts: 8,740	So now you're adding y twice in the RK function itself. You should only add it once (which is what you had when you started). Also F should return -y no matter what x is -- right now you're returning "?" most of the time from F, and who knows what's going to happen when you start adding "?" to things.
tabstop is offline

11-13-2008, 11:00 AM	#11
nickbrown05 Registered User Join Date: Nov 2008 Posts: 6	tabstop EST KING KING KING.............................................. ......................................... IT EFFING WORKS, for the first time I can see a bond between man and machine no longer will I live in fear of being snack attacked by the terminator
nickbrown05 is offline

Thread Tools
Show Printable Version Email this Page
Display Modes
Linear Mode Switch to Hybrid Mode Switch to Threaded Mode

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Need help understanding a problem	dnguyen1022	C++ Programming	2	04-29-2009 04:21 PM
Memory problem with Borland C 3.1	AZ1699	C Programming	16	11-16-2007 11:22 AM
Someone having same problem with Code Block?	ofayto	C++ Programming	1	07-12-2007 08:38 AM
A question related to strcmp	meili100	C++ Programming	6	07-07-2007 02:51 PM
WS_POPUP, continuation of old problem	blurrymadness	Windows Programming	1	04-20-2007 06:54 PM