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

---

Quote:

---

Originally Posted by civil25

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.

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.

> 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"

Quote:

---

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.

Quote:

---

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.

Quote:

---

Originally Posted by civil25

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

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! :)