printf() fixing pointer problems

Question

I'm trying to write a numerical integrator using a series expansion, however it currently fails on the second step forward in time. This gets fixed if I use printf to print out any of the variables, but I can't find the reason. I've already read many posts about similar scenarios, but none seem to answer my problem. If I print out the value of *t0, *Theta0 and *Theta1 within the function Integrate_One_Step(...), they all have the correct values and show t0, Theta0 and Theta1 to update with each step. Can someone tell me what I'm missing? See the code below. Thank you.

void main()
{
void Integrate_One_Step(double *Theta0, double *Theta1, double *t0, double t_final, double a, double epsilon, double gamma);

int i;
double a,epsilon,gamma,t0,t_final,Theta0,Theta1, tOld;

clock_t tic = clock();

epsilon = 0.1; a = 0.5; gamma = 0.07; t0 = 1e0; tOld = t0; t_final = 10;

Theta0 = 1.753317649649940; Theta1 = 1.759074746790801;

do{
    //printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
    //printf("tstart = %18.15f\n", t0);
    Integrate_One_Step(&Theta0,&Theta1,&t0,t_final,a,epsilon,gamma);
    //printf("tfinish = %18.15f\n", t0);
    //printf("theta = %18.15f\n", Theta0);
    //printf("tOld = %18.15f\n", tOld);
    //printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
    if(t0 <= tOld){
        printf("\nTime Step became zero\n");
        printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
        break;
    }
    else { tOld = t0; }
    if(Theta0 != Theta0){
        printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
        break;
    }
    if(Theta1 != Theta1){
        printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
        break;
    }
} while(t0 < t_final);

clock_t toc = clock();
printf("Elapsed: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);
}

Edit: I've posted the full code below for anyone who wishes to try and run it themselves. Thanks.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define OrderArray 25
#define HornerOrder 30

    void SeriesMult(double A[], double B[], double C[], int length);
    void Array_Division_By_Constant(double A[],double B[], double constant, int length);
    void Array_Multiply_By_Constant(double A[],double B[], double constant, int length);
    void ArrayAssign(double A[], double B[], int length);
    void Series_Add(double A[],double B[], double C[], int length);

int Time_Steps_Used[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double Fac[OrderArray],SinFac[HornerOrder],CosFac[HornerOrder],time_step[] = {0.005,0.01,0.03,0.06,0.1,0.2,0.4,0.6,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7};

void main()
{
    void Integrate_One_Step(double *Theta0, double *Theta1, double *t0, double t_final, double a, double epsilon, double gamma);

    int i;
    double a,epsilon,gamma,t0,t_final,Theta0,Theta1, tOld;

    clock_t tic = clock();

    epsilon = 0.1; a = 0.5; gamma = 0.07; t0 = 1e0; tOld = t0; t_final = 10;

    //Defines Fac[], SinFac[] and CosFac[] which are arrays of factorials
    Fac[0] = 1e0; SinFac[0] = 1;
    for (i = 1;i<= OrderArray-1;i++)    Fac[i] = i*Fac[i-1];
    for(i=1;i<=HornerOrder-1;i++){ SinFac[i] = -(2*i)*((2*i)+1); }
    for(i=0;i<=HornerOrder-1;i++){ CosFac[i] = -(2*(i+1))*((2*i)+1); }

    Theta0 = 1.753317649649940; Theta1 = 1.759074746790801;

    do{
        //printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
        //printf("tstart = %18.15f\n", t0);
        Integrate_One_Step(&Theta0,&Theta1,&t0,t_final,a,epsilon,gamma);
        //printf("tfinish = %18.15f\n", t0);
        //printf("theta = %18.15f\n", Theta0);
        //printf("tOld = %18.15f\n", tOld);
        //printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
        if(t0 <= tOld){
            printf("\nTime Step became zero\n");
            printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
            break;
        }else{ tOld = t0; }
        if(Theta0 != Theta0){
            printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
            break;
        }
        if(Theta1 != Theta1){
            printf("t = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);
            break;
        }
        //printf("tOld = %18.15f\n", tOld);
    } while(t0 < t_final);

    printf("\n\nt = %18.15f Theta = %18.15f dot_Theta = %18.15f\n", t0, Theta0,Theta1);

    printf("Time steps used: ");
    for(i = 0;i<=17;i++)    printf("%i ",Time_Steps_Used[i]);
    printf("\n\n");

    clock_t toc = clock();
    printf("Elapsed: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);
}

void Integrate_One_Step(double *Theta0, double *Theta1, double *t0, double t_final, double a, double epsilon, double gamma)
{
    double FGH_Nth_Position_Calc(double a, double epsilon, double gamma, double t0, double theta0, double X[], int Position, double dotTheta[], double ddotTheta[]);

    int i,k;
    double Theta[OrderArray],X[OrderArray],dotTheta[OrderArray-1],ddotTheta[OrderArray-2],/*FGH[OrderArray-2],*/    Thet,dotThet,ddotThet,EqnValue,Tol = pow(10,-12), t_inc = 0;

    //printf("At the start of step t = %18.15f \n",*t0);
    Theta[0] = *Theta0; Theta[1] = *Theta1; dotTheta[0] = Theta[1]; ddotTheta[0] = 0;
    X[0] = 0; //X is Theta without the constant term
    for(i=1;i<=OrderArray-1;i++){ X[i] = Theta[i];}

    for(i=0;i<=OrderArray-3;i++){
        //Calculate Theta i+2
        Theta[i+2] = FGH_Nth_Position_Calc(a,epsilon,gamma,(*t0),Theta[0],X,i,dotTheta,ddotTheta);
        //update X
        X[i+2] = Theta[i+2];
        //Update Theta' and Theta''
        dotTheta[i+1] = (i+2)*Theta[i+2];
        //ddotTheta[i] = (i+1)*(i+2)*Theta[i+2]; Slightly more multiplication, unnecessary.
        ddotTheta[i] = (i+1)*dotTheta[i+1];
    }
    //printf("Theta series: ");
    //for(i = 0;i<=OrderArray-1;i++)    printf("%18.15f ",Theta[i]);
    //printf("\n\n");

    //Calculate dotTheta and ddotTheta, then check if the step works
    for(k=0;k<=17;k++){
        ddotThet = 0.0; dotThet = 0.0; Thet = Theta[OrderArray-1];
        if(time_step[k] < (t_final - *t0)){
            for(i = OrderArray-2;i>=0;i--){
                ddotThet = ddotThet*time_step[k] + dotThet; dotThet = dotThet*time_step[k] + Thet; Thet = Thet*time_step[k] + Theta[i];
            }
            ddotThet *= 2e0;
            EqnValue = ddotThet + (((epsilon*sin(*t0+time_step[k]))/(1 + epsilon*cos(*t0+time_step[k])))+gamma)*dotThet + (a/(1 + epsilon*cos(*t0+time_step[k])))*sin(Thet);
            if(fabs(EqnValue) > Tol){ break; 
            }else{ t_inc = time_step[k]; *Theta0 = Thet; *Theta1 = dotThet; if(k>0){Time_Steps_Used[k-1] -= 1;} Time_Steps_Used[k] += 1; }
        }else{
            for(i = OrderArray-2;i>=0;i--){
                ddotThet = ddotThet*(t_final-*t0) + dotThet; dotThet = dotThet*(t_final-*t0) + Thet; Thet = Thet*(t_final-*t0) + Theta[i];
            }
            ddotThet *= 2e0;
            EqnValue = ddotThet + (((epsilon*sin(t_final))/(1 + epsilon*cos(t_final)))+gamma)*dotThet + (a/(1 + epsilon*cos(t_final)))*sin(Thet);
            if(fabs(EqnValue) > Tol){ break; }else{ t_inc = (t_final - *t0); *Theta0 = Thet; *Theta1 = dotThet; }
        }
    }
    //printf("After equation satisfied t0 = %18.15f\n",*t0);
    //printf("t_inc = %18.15f \n",t_inc);
    printf("theta = %18.15f dot_theta = %18.15f\n\n",*Theta0,*Theta1);
    //Add time step to t0
    *t0 = *t0 + t_inc;
    //printf("After incremented t = %18.15f\n",*t0);
}

double FGH_Nth_Position_Calc(double a, double epsilon, double gamma, double t0, double theta0, double X[], int Position, double dotTheta[], double ddotTheta[])
{
    void CosT_Calc(double coss[], double t0, int length);
    void FCalc(double epsilon, double Cost[], double ddotTheta[], double F[], int length);
    void HCalc(double a, double X[], double H[], double theta0, int Position);
    void GCalc(double epsilon, double gamma, double t0, double Cost[], double G[], double dotTheta[], int Position);

    int i;
    double F[Position+1],H[Position+1],G[Position+1],Cost[Position+1],FGH[Position+1],Divider;

    //for(i = 0;i<=Position;i++)    printf("%18.15f ",dotTheta[i]);
    //printf("\n\n");
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",ddotTheta[i]);
    //printf("\n\n");
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",ddotTheta2[i]);
    //printf("\n\n");

    //Calculates coefficicients of cos(t) where t is the total time after the step is taken.
    CosT_Calc(Cost,t0,Position+1);
    //for(i=0;i<=Position;i++)  printf("%18.15f ",Cost[i]); printf("\n");

    //Calculates F = epsilon*cos(t - t0)* \theta''
    FCalc(epsilon,Cost, ddotTheta, F, Position);
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",F[i]);
    //printf("\n\n");
    //Calculates G = [epsilon*sin(t) + gamma*(1+epsilon*cos(t)) ]*theta'
    GCalc(epsilon,gamma,t0,Cost,G,dotTheta,Position);
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",G[i]);
    //printf("\n\n");
    //Calculates H = a*sin(Theta)
    HCalc(a,X, H,theta0,Position);
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",H[i]);
    //printf("\n\n");

    Series_Add(F,G,FGH,Position+1);
    Series_Add(FGH,H,FGH,Position+1);
    Divider = (Position+1)*(Position+2)*(1+epsilon*cos(t0)); //printf("Divider = %18.15f\n",Divider);
    Array_Division_By_Constant(FGH,FGH,Divider, Position+1);
    return(-FGH[Position]);
}

void GCalc(double epsilon, double gamma, double t0, double Cost[], double G[], double dotTheta[], int Position)
{
    void SinT_Calc(double Sinn[], double t0, int length);

    int i;
    double Sint[Position+1],epsilSint[Position+1],epsilCostPlusOne[Position+1],g[Position+1];

    //Calculates sin(t)
    SinT_Calc(Sint,t0,Position+1);
    //printf("Sin(t) series");
    //for(i=0;i<=Position;i++)  printf("%18.15f ",Sint[i]); printf("\n");
    //multiplies sin(t) and cos(t) by epsilon
    Array_Multiply_By_Constant(epsilSint,Sint,epsilon,Position+1);
    Array_Multiply_By_Constant(epsilCostPlusOne,Cost,epsilon,Position+1);
    //Calculates gamma*(1 + epsilon*cos(t))
    epsilCostPlusOne[0] = epsilCostPlusOne[0]+1;
    Array_Multiply_By_Constant(epsilCostPlusOne,epsilCostPlusOne,gamma,Position+1);
    //Calculates g = [epsilon*sin(t) + gamma*(1+epsilon*cos(t)) ]
    Series_Add(epsilSint,epsilCostPlusOne,g,Position+1); //This appears to be correct when comparing with maple
    //for(i = 0;i<=Position;i++)    printf("%18.15f ",g[i]);
    //printf("\n\n");
    //Multiples g by theta'
    SeriesMult(g,dotTheta,G,Position+1);

}

void HCalc(double a, double X[], double H[], double theta0, int Position)
{
    void SinTheta_Calc(double CosX[], double SinX[], double SinTheta[], double theta0, int length);
    void CosX_Calc(double Xsquared[], double CosTemp[], int length);
    void SinX_Calc(double X[], double Xsquared[], double Sinn[], int length);
    double NthTerm(double A[], double B[], int n, int sizeA, int sizeB);

    int i;
    double Xsquared[Position+1],CosX[Position+1],SinX[Position+1],SinTheta[Position+1];
    //Calculates Xsquared
    for(i=0;i<=Position;i++){
        Xsquared[i] = NthTerm(X,X,i,Position,Position);
    }
    //Calulate Cos(X), Sin(X)
    CosX_Calc(Xsquared,CosX,Position+1); SinX_Calc(X,Xsquared,SinX,Position+1);
    //Calculate Sin(Theta) from Sin(X) and Cos(X)
    SinTheta_Calc(CosX,SinX,SinTheta,theta0, Position+1);
    //H = a*SinTheta
    Array_Multiply_By_Constant(H,SinTheta, a, Position+1);

}

void SinTheta_Calc(double CosX[], double SinX[], double SinTheta[], double theta0, int length)
{
    int i;
    double CosXTemp[length],SinXTemp[length], S_Theta0, C_Theta0;
    S_Theta0 = sin(theta0); C_Theta0 = cos(theta0);
    Array_Multiply_By_Constant(SinXTemp,SinX,C_Theta0,length);
    Array_Multiply_By_Constant(CosXTemp,CosX,S_Theta0,length);
    Series_Add(CosXTemp,SinXTemp,SinTheta,length);
}

void CosX_Calc(double Xsquared[], double CosTemp[], int length)
{
    //Calculates Cos(X) where X is the series for Theta, without the constant term, outputs the final series for Cos(X) as CosTemp
    int i,j,SubHornerOrder;
    double Xsquared_divided[length], Coss[length];

    if(length == 0){SubHornerOrder = 1; }
    else if(length%2 == 0) {SubHornerOrder = length/2; }
    else {SubHornerOrder = (length+1)/2; }

    Array_Division_By_Constant(CosTemp,Xsquared,CosFac[SubHornerOrder-1],length);
    CosTemp[0] = 1;
    for(i=(SubHornerOrder-1);i>=1;i--){
        Array_Division_By_Constant(Xsquared_divided,Xsquared,CosFac[i-1],length);
        SeriesMult(Xsquared_divided, CosTemp, Coss, length);
        ArrayAssign(CosTemp, Coss, length);
        CosTemp[0] = 1;
    }
}



void SinX_Calc(double X[], double Xsquared[], double Sinn[], int length)
{
    //Calculate Sin(X) where X is the series for Theta, without the constant term, outputs the final series for Sin(X) as Sinn
    int i,j, SubHornerOrder;
    double Xsquared_divided[length], SinTemp[length];

    if(length == 0){SubHornerOrder = 1; }
    else if(length%2 == 0) {SubHornerOrder = length/2; }
    else {SubHornerOrder = (length+1)/2; }

    Array_Division_By_Constant(SinTemp,Xsquared,SinFac[SubHornerOrder-1],length);
    SinTemp[0] = 1;
    for(i=(SubHornerOrder-1);i>=2;i--){
        Array_Division_By_Constant(Xsquared_divided,Xsquared,SinFac[i-1],length);
        SeriesMult(Xsquared_divided, SinTemp, Sinn, length);
        ArrayAssign(SinTemp, Sinn, length);
        SinTemp[0] = 1;
    }
    SeriesMult(X,SinTemp,Sinn,length);
}

void FCalc(double epsilon, double Cost[], double ddotTheta[], double F[], int Position)
{
    double epsilCost[Position+1];
    //multiplies Cos(t) by epsilon
    Array_Multiply_By_Constant(epsilCost,Cost,epsilon,Position+1);
    //Calculates F = epsilon*cos(t - t0)* \theta''
    SeriesMult(epsilCost, ddotTheta, F, Position+1);
}

void SinT_Calc(double sinn[], double t0, int length)
{ // Sets up Taylor series coefficients for sin(t) = cos(t0)sin(t-t0) + sin(t0)cos(t-t0)
    int i;
    double  st0, ct0;
    int m1n(int i);

    st0 = sin(t0); ct0 = cos(t0);

    for (i = 0; i <= length; i++)
    {
        if (i%2 == 0)
            sinn[i] = st0*m1n(i/2)/Fac[i];
        else
            sinn[i] = ct0*m1n((i-1)/2)/Fac[i];
    }
}

void CosT_Calc(double coss[], double t0, int length)
{ // Sets up Taylor series coefficients for cos(t) = cos(t0)cos(t-t0) - sin(t0)sin(t-t0)
    int i;
    double  st0, ct0;
    int m1n(int i);

    st0 = sin(t0); ct0 = cos(t0);

    for (i = 0; i <= length; i++)
    {
        if (i%2 == 0)
            coss[i] = ct0*m1n(i/2)/Fac[i];
        else
            coss[i] = -st0*m1n((i-1)/2)/Fac[i];
    }
}

int m1n(int i)
{
    if (i%2 == 0)   return(1);
    else        return(-1);
}

double NthTerm(double A[], double B[], int n, int sizeA, int sizeB)
{
    //Multiplies returns the n-th term when 2 arrays are multiplied
    double nth_term = 0;
    int j;
    for(j=0; j <= n;j++){
        if (j < sizeA){
            if ( (n-j)< sizeB){
                nth_term += A[j]*B[n-j];
            }
        }
    }
    return(nth_term);
}

void Series_Add(double A[],double B[], double C[], int length)
{
    int i;
    for(i=0;i<=length-1;i++){
        C[i] = A[i] + B[i];
    }
}

void SeriesMult(double A[], double B[], double C[], int length)
{
    //Multiplies two arrays, returns nothing but stores the answer in the array C
    int i, j;
    double c;
    for(i=0; i<=length-1; i++){
        c = 0;
        for(j=0; j<=i; j++){
            c = c + A[j]*B[i-j];
        }
        C[i] = c;
    }
}

void Array_Multiply_By_Constant(double A[],double B[], double constant, int length)
{
    int i;
    for(i=0;i<=length-1;i++){
        A[i] = B[i]*constant;
    }
}

void Array_Division_By_Constant(double A[],double B[], double constant, int length)
{
    int i;
    for(i=0;i<=length-1;i++){
        A[i] = B[i]/constant;
    }
}

void ArrayAssign(double A[], double B[], int length)
{
    int i;
    for(i=0;i<=length-1;i++){
        A[i] = B[i];
    }
}

Answer 1

You have buffer overflow errors in places such as:

void SinT_Calc(double sinn[], double t0, int length)
{ // Sets up Taylor series coefficients for sin(t) = cos(t0)sin(t-t0) + sin(t0)cos(t-t0)
    int i;
    double  st0, ct0;
    int m1n(int i);

    st0 = sin(t0); ct0 = cos(t0);

    for (i = 0; i <= length; i++) // <=== BUFFER OVERFLOW, should be strict <
    {
        if (i%2 == 0)
            sinn[i] = st0*m1n(i/2)/Fac[i];
        else
            sinn[i] = ct0*m1n((i-1)/2)/Fac[i];
    }
}

If I had one other unrelated piece of advice, it would be to never use dynamically sized arrays on the stack (actually, it's surprising this even compiles):

double Xsquared[Position+1];

Instead use

double* Xsquared = (double*)malloc((Position + 1) * sizeof(double));
//...
free(Xsquared);