printf() проблемы с фиксацией указателя

Question

Я пытаюсь написать числовой интегратор с использованием расширения серии, однако в настоящее время он не работает на втором шаге вперед во времени. Это фиксируется, если я использую printf для печати любой из переменных, но я не могу найти причину. Я уже прочитал много сообщений о подобных сценариях, но никто, кажется, не ответил на мою проблему. Если я распечатаю значение * t0, * Theta0 и * Theta1 в функции Integrate_One_Step (...), все они имеют правильные значения и показывают, что t0, Theta0 и Theta1 обновляются с каждым шагом. Может ли кто-нибудь сказать мне, что мне не хватает? См. Код ниже. Спасибо.

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: Я отправил полный код ниже для тех, кто хочет попробовать и запустить его самостоятельно. Спасибо.

#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

У вас есть буфер ошибок переполнения в таких местах, как:

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

Если бы я был еще один несвязанный совет, было бы никогда не использовать динамически размера массивов в стеке (на самом деле, это удивительно это даже компилируется):

double Xsquared[Position+1]; 

Вместо этого используйте

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