Runge Kutta IV project problem

[achikha]

Hey guys, i had a project to do, , this one I figured code but:1. It's not working properly
2. Mechanic energy isn't constant(it should be )

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define g 9.81
#define PI 4.*atan(1.)
#define MAXN 2


double t0, alfa0, omega0, R, t, NS;




    void vrk4( double x0, double y0[], double h, int n, void (*fun)(double, double*, double*), double y1[] )
{
	int		i;
	double	k1[MAXN], k2[MAXN], k3[MAXN], k4[MAXN];
	double	ytmp[MAXN];
	fun( x0, y0, k1);
	for ( i=0; i<n; ++i)
	{
		k1[i] *= h;
		ytmp[i] = y0[i] + k1[i]/2.0;
	}


	fun( x0+h/2.0, ytmp, k2);
	for ( i=0; i<n; ++i)
	{
		k2[i] *= h;
		ytmp[i] = y0[i] + k2[i]/2.0;
	}


	fun( x0+h/2.0, ytmp, k3);
	for ( i=0; i<n; ++i)
	{
		k3[i] *= h;
		ytmp[i] = y0[i] + k3[i];
	}


	fun( x0+h, ytmp, k4);
	for ( i=0; i<n; ++i)
		k4[i] *= h;


	for ( i=0; i<n; ++i)
		y1[i] = y0[i] + (k1[i] + 2.*k2[i] + 2.*k3[i] + k4[i])/6.;


}


double fun(double x0, double *y0, double *F)
{
    F[0]=y0[1];
    F[1]=g*sin(y0[0])/R;
}


int main()
{
    double m=1;
    R=0.5;
    alfa0=0.;
    omega0=(1./3.)*PI*R;
    t=0.1;
    NS=100;
    double y0[2]={alfa0,omega0};
    double x0=0.;
    double y1[2];
    double ss;
    int p,j;
    FILE *f;
    f=fopen("Results.txt", "wt");
    fprintf(f,"t\talfa\tomega\n");
    fprintf(f,"%lg\t%lg\t%lg\t%lg\n", x0,y0[0],y0[1],m*pow(y0[1]*R,2.)/2.+m*g*R*(1.-cos(y0[0])));
    printf("t\talfa\tomega\n");
    for(p=0;p<NS;p++)
    {
        ss=1.-cos(y1[0]);
        vrk4(x0,y0,t,2,fun,y1);
y0[0]=y1[0];
y0[1]=y1[1];


       // for(j=0;j<2;j++)
       // {
       //     y0[j]=y1[j];
       // }
        x0+=t;
        printf("%lg\t%lg\t%lg\n",x0,y1[0],y1[1]);
        fprintf(f,"%lg\t%lg\t%lg\t%lg\n", x0,y1[0],y1[1],m*pow(y1[1]*R,2.)/2.+m*g*R*(ss));
    }
    fclose(f);
}

Can anybody help me please?

[oogabooga]

You probably need to give more information to get some help.
"It's not working properly" isn't enough.

[achikha]

You probably need to give more information to get some help.
"It's not working properly" isn't enough.

okay, You may be right, sorry.
for what i think it's the for loop in the main function that's miswritten.
The program finds the first derivative after t, but fails to find second one. I don't know where exactly i've made a mistake but seems that in the for loop in main().

[vart]

The program finds the first derivative after t, but fails to find second one. I don't know where exactly i've made a mistake but seems that in the for loop in main().

So start writing unit testing.

Create simple main that just uses you function to create a first derivative and print the results.

When this is working fine.

Add the second call - which will calculate the second derivative... using the previous output as input.

Check teh results again...

Only after this 2 steps are working correctly - continue with the actual task.

[achikha]

yes, I still get -1.#IND00 in the 2nd derivative, having changed the code into :

#include <stdio.h>#include <stdlib.h>
#include <math.h>
#define g 9.81
#define PI 4.*atan(1.)
#define MAXN 2




double t0, alfa0, omega0, R, t, NS;








    void vrk4( double x0, double y0[], double h, int n, void (*fun)(double, double*, double*), double y1[] )
{
    int     i;
    double  k1[MAXN], k2[MAXN], k3[MAXN], k4[MAXN];
    double  ytmp[MAXN];
    fun( x0, y0, k1);
    for ( i=0; i<n; ++i)
    {
        k1[i] *= h;
        ytmp[i] = y0[i] + k1[i]/2.0;
    }




    fun( x0+h/2.0, ytmp, k2);
    for ( i=0; i<n; ++i)
    {
        k2[i] *= h;
        ytmp[i] = y0[i] + k2[i]/2.0;
    }




    fun( x0+h/2.0, ytmp, k3);
    for ( i=0; i<n; ++i)
    {
        k3[i] *= h;
        ytmp[i] = y0[i] + k3[i];
    }




    fun( x0+h, ytmp, k4);
    for ( i=0; i<n; ++i)
        k4[i] *= h;




    for ( i=0; i<n; ++i)
        y1[i] = y0[i] + (k1[i] + 2.*k2[i] + 2.*k3[i] + k4[i])/6.;




}




double fun(double x0, double *y0, double *F)
{


    F[0]=y0[1];
    F[1]=g*sin(y0[0])/R;


}




int main()
{
    double m=1;
    R=0.5;
    alfa0=2.;
    omega0=(1./9.)*PI*R;
    t=0.1;
    NS=100;
    double F[2];
    double y0[2]={alfa0,omega0};
    double x0=0.;
    double y1[2];
    double ss;
    int p,j;
    FILE *f;
    f=fopen("Results.txt", "wt");
    fprintf(f,"t\talfa\tomega\n");
    fprintf(f,"%lg\t%lg\t%lg\t%lg\n", x0,y0[0],y0[1],m*pow(y0[1]*R,2.)/2.+m*g*R*(1.-cos(y0[0])));
    printf("t\talfa\tomega\n");
    printf("%lf\t\t%lf\t\t%lf\n",x0,y0[0],y0[1]);
    x0+=t;
    vrk4(x0,y0,t,2,fun,y1);
    printf("%lf\t\t%lf\t\t%lf\n",x0,y1[0],y1[1]);
    y0[0]=y1[0];
    y0[1]=y1[1];
    x0+=t;
    vrk4(x0,y0,t,2,fun,y1);
    printf("%lf\t\t%lf\t\t%lf\n",x0,y1[0],y1[1]);




//    for(p=0;p<NS;p++)
//    {
//        ss=1.-cos(y1[0]);
//        vrk4(x0,y0,t,2,fun,y1);
//y0[0]=y1[0];
//y0[1]=y1[1];
//
//
//       // for(j=0;j<2;j++)
//       // {
//       //     y0[j]=y1[j];
//       // }
//        x0+=t;
//        printf("%lg\t%lg\t%lg\n",x0,y1[0],y1[1]);
//        fprintf(f,"%lg\t%lg\t%lg\t%lg\n", x0,y1[0],y1[1],m*pow(y1[1]*R,2.)/2.+m*g*R*(ss));
//    }
    fclose(f);
}

[rcgldr]

It's difficult to read the image. I would think that you're supposed to show position(t), velocity(t), acceleration(t), energy(t), but I'm not sure what alpha, omega, and the third value you're supposed to calculate is.

You're rk4 function seems off a bit, k1 should not be divided by 2 (it's an euler step), while both k2 and k3 should be divide by 2 (they're the mid-point steps). k4 seems ok (it's not divided by 2). update - ok, I see that ytmp is being used for the next k#, in that case it's OK.

This pdf file on using rk4 for orbits may help, as it shows the basic logic for dealing with position, velocity, and acceleration.

http://spiff.rit.edu/richmond/nbody/...ungeKutta4.pdf

[achikha]

i'm pretty sure about the function rk4 , it's a vector rk4 (we got it from our teacher), and the values to calculate is angular velocity (t), angle (t), angular velocity(angle)

[rcgldr]

i'm pretty sure about the function rk4 , it's a vector rk4 (we got it from our teacher), and the values to calculate is angular velocity (t), angle (t), angular velocity(angle)

I updated my previous post, I didn't notice the order of your operations with ytmp. It might help if you also output angular position(t) = theta.

I don't understand fun(), you set y[0] to angular acceleration in main, but then in fun(), you calculate F[1] to g sin(y[0]) / R.

Shouldn't angular acceleration be a function of theta?

Assuming the conventional direction for positive angular velocity ω is counter clockwise, then R pointed directly left corresponds to θ = -π / 2, R pointed directly down corresponds to θ = 0, and R pointed directly right corresponds to θ = π / 2. In this case, it seems that angular acceleration should be α = -g sin(θ) / R?

[achikha]

y_tmp is just to get the derivative after time step y1[0] and y1[1] so alfa and omega
differential eqns (initial value problems) are: alfa'=omega, alfa(0)=alfa0, omega'=g*sin(alfa)/R, omega(0)=omega0;

[achikha]

hmm so y0[0] is alfa0 (starting position) and y0[1] is starting angular velocity, then y1[0] final pos y1[1] final velocity
edit: in here there is no need to use angular acceleration , the problem is somewhere in the fun function (dividing 0/0 i don't know how) or in the for loop in main :/

[achikha]

okay, that makes sense, but how to put it in code :P?

[vart]

the problem is somewhere in the fun function (dividing 0/0 i don't know how) :/

so add trace output in the fun function and see all the incoming arguments before calculations...

or use the debugger

[rcgldr]

y_tmp is just to get the derivative after time step y1[0] and y1[1] so alfa and omega
differential eqns (initial value problems) are: alfa'=omega, alfa(0)=alfa0, omega'=g*sin(alfa)/R, omega(0)=omega0;

You're not using the normal terminology for angular stuff, which is confusing to me.

angular position = theta = θ
angular velocity = omega = ω
angular acceleration = alpha = α

[achikha]

sorry for that, could be theta :P

[achikha]

can't somebody compile it and check why i get -1.#IND00 pleaaase?