4 Replies - 1081 Views - Last Post: 08 April 2011 - 01:33 AM Rate Topic: -----

#1 nzm  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 3
  • Joined: 03-April 11

the output of this program is SOL[1][0]=-1.#IND000000e+000.

Posted 03 April 2011 - 07:40 PM

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <usclkc.h>
#define true 1
#define false 0

/* GLOBAL VARIABLES DECLARATION*/

int N,TEST,MAXD,TESTS,STIFF;
double X,XEND,TOL,H;
double Y[9][9];
int D[9],J[9];

/*new add for global declaration*/

double HMIN,REL,AB,ED[9],RD[9],P[9][9],PSI[21],F[9],MAX,DELTA,YP[9],MAXE,XMAXE;
double F1[9][14],G[15][17],ER[9],JACBN[9][9],B[9],ERS[9],HH;
int FIRSTSTEP,STORE[9],FR,KM1,RR[9],K[9],S[9],NEWJAC,JUMP,CNT[9];
int ST[9],MJ[9],NS,MKP1,T,TT,JMAX,PS[9],R1; 

/**********************************************************************/
/*                         STARTS MAIN PROGRAM                        */
/**********************************************************************/

main()

{

	void INTEGRATE2();
	TOL=1.0E-7; /* TOL=tolerance limit */

	N=1;/* number of equations*/
	TEST=1;
	TOL=TOL*0.1;/* tolerance limit*/
	Y[1][0]=(double)15.e0;/* initial condition for test problem*/
	
	D[1]=1;/* order of the test problem*/
	

	X=(double)2.e0; /* initial value */
	XEND=(double)(25.e0);/* end point*/
	
	H=(double)1.e-4;
	J[1]=0;
	
    MAXD=1;

	INTEGRATE2(); /* Enter subroutine Integrate2*/
}

/*********************************************************************/
/*                     END OF MAIN PROGRAM                           */
/*********************************************************************/

/*********************************************************************/
/*                  STARTS SUBROUTINE INTEGRATE2                     */
/*********************************************************************/

void INTEGRATE2()

{

/* FUNCTIONS DECLARATION*/

void FCN(),SOLN(),STSIZE(),EVALUATE(),DERV();


/* VARIABLES DECLARATION*/

double TEMP,TEMP1,TEMP2,HOLD,XOLD,ERRMAX,KT,EPS,
       RAT,SIGMA,SUM,E1,E2,TEMP3,TEMP4;

double E[16],MAXEG[3],DC[15][17],KH[5],ABSE[9],F2[9][4];
  
int    I,L,MAXJ,IM1,KP1,IP1,IND,KOLD,ERRORS,STEPS,
       K2,SS,GMAX,JP1,JP2,FACT,EVALE2;
 
int    SUCSTP,FA,TOTALSTP,BOOL[9],DUMMY;

double N10,N11,N12,N13,DEL;

 FILE*f1;
	f1=fopen("C:PROFPRGRM.txt","w");

	if (TEST==2) REL=(double)0.e0; 
	else 
	    REL=(double)1.e0;

	if (TEST==3) AB=(double)0.e0; 
	else 
	    AB=(double)1.e0;
	
	if (TOL>1.0e-9) 
			FR=3; 
		else 
            FR=4;

        MAXE=(double)0.e0;
	    XMAXE=(double)0.e0;

	    HMIN = HOLD=(double)0.e0;
	    NS=0;

	 // TESTS=0; 
	    JMAX=0; 
	    EPS=1.0e-48; 
	    JUMP=false;
	    
		for (I=1;I<=N;I++)
		{
			BOOL[I]=false; 
	    	ERS[I]=(double)0.e0;
		}

	    for (I=1;I<=4;I++)
		    ER[I]=0.e0;
			ED[I]=0.e0;
		    
	    FIRSTSTEP=true;
	    STEPS=0;
		SUCSTP=0;
		FA=0;
		TOTALSTP=0;
	   	  
	    EVALUATE(1,N,Y,F,2);
	
L45: 		for (I=1;I<=N;I++) 
			Y[I][D[I]]=F[I];

L43: 		for (I=1;I<=N;I++)
			{
    			MJ[I]=D[I]-J[I];
     		    

		    if (MJ[I]==0) 
				S[I]=0;
			else if (MJ[I]==1) 
				S[I]=1;
			else if (MJ[I] <= 4) 
				S[I]=2; 
			else 
				S[I]=3;

		        RR[I]=J[I];

			}/*end I*/

	       if(H==0)STSIZE(J,D,Y,&H);
	       fprintf(f1,"Initial Stepsize=%10.5e\n\n",(double)H);
  
L4:     MKP1=0;
		for (I=1;I<=N;I++)
		{
			K[I]=J[I]+1;

			if (K[I]>MKP1) MKP1=K[I];
	    }/*end I*/

	    MKP1=MKP1+1;
    	DERV(F1,Y,MJ,K);
	    
L7:    	for(I=1;I<=MKP1;I++)
	
		PSI[I]=(double)0.e0;
		
/*********************************************************************/
/*                      COMPUTE CHANGE COEFFICIENTS                  */
/*********************************************************************/
L10:    if (H != HOLD) NS=0;
		
		if (NS >= MKP1) goto L11;
	
		if (NS<MKP1) NS=NS+1;

		if(MKP1>NS)
			for (I=MKP1;I>=NS+1;I--)
			{
			    PSI[I]=PSI[I-1]+H;
			}
				PSI[NS]=H*(double)(NS);
						
				for (I=NS;I<=MKP1;I++)
				{
		
					G[I][0]=PSI[I]*G[I-1][0];
				}	
				for (K2=1;K2<=N;K2++)
			
					if (MJ[K2] != 0)
					{
						GMAX=0;
						for (I=1;I<=N;I++)
						{
					    	L=K[I]+MJ[I];
						
							if (GMAX<L) 
								GMAX=L;
						}/*end I*/

					if (NS==1)
						{
							for (T=1;T<=GMAX+1;T++)
							{
								G[1][T-1]=G[0][T-1]*H/(double)(T);
								G[0][T]=G[1][T-1];
							}
						}/*end NS==1*/

					if ((MKP1>KOLD) && (NS>1))
						{
							G[1][GMAX]=H*G[1][GMAX-1]/(double)(GMAX+1);
							G[0][GMAX+1]=G[1][GMAX];

							for (I=2;I<=NS-1;I++)
							{
								T=GMAX+1-I;
								IM1=I-1;
								G[I][T]=PSI[I]*G[IM1][T]-((double)(T)*G[IM1][T+1]);
							}/*end I*/
						}/*end MKP1>KOLD*/

					if (NS==1) 
							L=2; 
						else 
							L=NS;

						for (I=L;I<=MKP1;I++)
						{
							IM1=I-1; 
							MAXJ= GMAX+1-I;
						    for (T=1;T<=MAXJ;T++)
							{
								G[I][T]=PSI[I]*G[IM1][T]-T*G[IM1][T+1];
							}
						}/*end I*/
						// goto L07;
					}/*end TT<N*/


		

L11:			XOLD=X;
				X=X+H;
				fprintf(f1,"X=%15.10e,\tH=%15.10e\n\n",X,H);
				printf("X=%15.10e,\tH=%15.10e\n\n",X,H);
				HOLD=H;
				KOLD=MKP1;

L14:  			HMIN=(double)3.e0;
				TOTALSTP=TOTALSTP+1;

	  
/************************************************************************/
/*                                  PREDICTION                          */
/************************************************************************/

for (I=1;I<=N;I++)
			 {
				 fprintf(f1,"K[%d]=%d\n\n",I,K[I]);
				 printf("K[%d]=%d\n\n",I,K[I]);
				for (K2=0;K2<=MJ[I];K2++)
				{
				  SIGMA=0.0;
				  if (K2==0) SUM=(double)0.e0;
				  else 
				  SUM=Y[I][MJ[I]-1];
				   for (T=0;T<=K[I]-1;T++)
						SIGMA=SIGMA + G[T][K2]*F1[I][T];
						  	if (K2<2)  goto L49;
												
							for (L=(K2-1);L>=1;L--)
							
								SUM=SUM*(H/L)+Y[I][MJ[I]-K2+L-1];

L49:	P[I][MJ[I]-K2]=SUM + SIGMA;/*predicted value*/
			fprintf(f1,"P[%d][%d]=%15.10e\n",I,MJ[I]-K2,P[I][MJ[I]-K2]);
				}/* end K2*/
							
						
					}/* end I*/


EVALUATE(1,N,P,F,0);
	for (I=1;I<=N;I++)
	if (J[I]==0)
	{
		ED[I]=F[I]-P[I][D[I]];
		for (L=1;L<=D[I];L++)
			P[I][MJ[I]-L]=P[I][MJ[I]-L]+ED[I]*G[K[I]][L]/G[K[I]][0];
	}
	    fprintf(f1,"P[%d][%d]=%15.10e\n",I,MJ[I]-L,P[I][MJ[I]-L]);						
		if(TT>=1) EVALUATE(1,N,P,F,1); /*try yang ni */


/***********************************************************************/
/*                                 ERROR COMPUTATION                   */
/***********************************************************************/


	ERRMAX=(double)0.e0;
	
	for(I=1;I<=N;I++){
		KM1=K[I]-1;

		TEMP=AB+REL*fabs(P[I][MJ[I]-S[I]]);
	if (J[I]==0)
	{
	E[I]=fabs((S[I])*G[KM1][S[I]+1]*ED[I]/G[K[I]][0])/TEMP;
	}
L08:	if (ERRMAX<E[I]){
			ERRMAX=E[I];
		
	}/*end I*/

	if (ERRMAX<TOL) 
	{
		IND=1;
	}
	else{
		  FA++;
		  IND=2;
		  goto L15;
	}

SUCSTP=SUCSTP+1; 
}
	
/***********************************************************************/
/*                                  CORRECTIONS                        */
/***********************************************************************/

	 EVALUATE(1,N,P,F,0);
/***************************************************************/
/*                    DIVIDED DIFFERENCES                      */
/***************************************************************/

L15:	if (IND==1){
		for (I=1;I<=N;I++){
	
			Y[I][D[I]]=F[I];

			for (L=1;L<=D[I];L++) 
				Y[I][D[I]-L]=P[I][D[I]-L];
			
			if (J[I]==0) 
				TEMP=F[I];
			else 
				TEMP=P[I][MJ[I]];

			if (FIRSTSTEP) 
				T=K[I]-1;
			else 
				T=K[I];

			for (L=0;L<=T;L++){
				TEMP1=(TEMP-F1[I][L])/PSI[L+1];
				F1[I][L]=TEMP; 
				TEMP=TEMP1;
			}/*end for statement L*/
			
			F1[I][T+1]=TEMP;
		}/*end I*/
		SOLN();
	}/*end IND==1*/
	else{	
		for (I=1;I<=N;I++){
			F2[I][3]=ED[I]/(G[K[I]-1][0]*PSI[K[I]]);

			if (K[I]>J[I]+1)
				F2[I][2]=F1[I][K[I]-1]+PSI[K[I]]*F2[I][3];

			if (K[I]>J[I]+2)
				F2[I][1]=F1[I][K[I]-2]+PSI[K[I]-1]*F2[I][2]; 
			
		}/*end IND==2*/
	}/*end I*/

if (FIRSTSTEP){
	        if (IND==1)
			{
		    	FIRSTSTEP=false;
			    for (I=1;I<=N;I++)
					ERRORS=0;
	           goto L27;
			}
				else 
		{
			if (ERRORS==0) 
				HMIN=(double)1.e0/(double)(1.e1);
			else if (ERRORS==1)
				HMIN=(double)1.e0/(double)(1.e3);
		}
				}/*end FIRSTSTEP*/

/**************************************************************************/
/*                              ORDER AND STEPSIZE                        */
/**************************************************************************/

MKP1=0;
		for (I=1;I<=N;I++)
	{	
		KP1=K[I]+1;
		JP1=J[I]+1;
		KM1=K[I]-1;
		JP2=JP1+1;

       if (K[I]>JP1+1)
	   	      SS=1;
	   else if (K[I]>JP1){
		   SS=2;}
	   else {
		   SS=3;}
   		for(L=SS;L<=3;L++)
		{
		T=KM1+L-2;
		if (IND==1)
			TEMP=F1[I][T];
			else 
			  TEMP=F2[I][L];
			if (J[I]==0)
				ER[L]=(double)fabs(S[I]*G[T-1][S[I]+1]*TEMP);
			
		}/*end L*/
		if ((K[I]==JP2) && (ER[2]<(0.5*ER[3])))
			 K[I]=KM1;
			
		if (K[I]>JP2)
			{
				if (ER[1]>ER[2])
					TEMP1=ER[1];
				else 
					TEMP1=ER[2];
				if (TEMP1<ER[3]) 
			    	K[I]=KM1;
			
			}
		if(J[I]==0)
			N13=K[I]<12;
			else 
			N13=K[I]<6;
		if ((! BOOL [I]) && (K[I] != KM1) && (IND==1) && (N13))
		{
			RAT=(double)(2.e0);
			K[I]=KP1;
			goto L12;
		} /*end BOOL[I]*/
		else 
			BOOL[I]=true;
		if (IND==1){
			if((NS >= KP1) && (K[I] != KM1) && (J[I]==0 && K[I]<12) /*|| (J[I]>0 && K[I]<6)*/) 
{
  if (J[I]==0) 
	ER[4]=(double)fabs(S[I]*G[K[I]][S[I]+1]*F1[I][KP1]);
	else 
    ER[4]=(double)fabs((DC[K[I]][J[I]]/DC[KP1][J[I]]*G[KP1][S[I]]-G[K[I]][S[I]])*PSI[KP1]*F1[I][KP1]);
		if (K[I]==JP1)
		{
		 if (ER[4]<(double)(0.5)*ER[3]) 
			K[I]=KP1;
		}
		else
		{
		  	if ((/*(J[I]>0 && ER[4]<ER[3]) || */ (J[I]==0 && ER[4]<ER[3])) && ER[3]<ER[2])
			K[I]=KP1;
		    else if (ER[2]<ER[3] && ER[2]<ER[4])
		        	K[I]=KM1;
		}/*end K[I]>1*/
	}/*end IND=1*/

}/*end IND=1*/

	TEMP2=AB+REL*(double)fabs(P[I][MJ[I]-S[I]]);
	T=K[I]-KM1+2;

	if (ER[T]<EPS) 
		RAT=(double)(4.e0);
	else 
		RAT=(double)pow(TOL*TEMP2/ER[T],(double)1.e0/(double)(K[I]+S[I]));
	

L12:	if (HMIN>RAT) HMIN=RAT;
		if (MKP1<K[I]) MKP1=K[I];

	}/*end I*/
	
	MKP1=MKP1+1;

L22:	if (IND==2){
		ERRORS=ERRORS+1;
	    
		if (! FIRSTSTEP){
			if (HMIN <=(double)0.5e0) 
				HMIN=0.4;
			else if (HMIN>(double)0.5e0 && HMIN <=(double)0.8e0) 
				HMIN=(double)0.5e0;
				else 
				HMIN=(double)0.7e0;
				}

		H=HMIN*HOLD;
		X=XOLD;

		if (NS==0) 
			T=2;
		else 
			T=NS+1;
	
		for(T=NS+1;T<=KOLD;T++) 
			PSI[T-1]=PSI[T]-HOLD;
	 

		if (ERRORS==12){
			printf("TOO MANY ERRORS %20.5e %d %d\n",X,0,5);
			goto g302;
		}
		else if (ERRORS>3){
			H=H/(double)(10);
			goto L4;
		}
	}/*end IND==2*/
	else{
		HMIN=HMIN*(double)0.8e0;
		
		if (HMIN >=(double)(2)) 
		
			H=(double)(2.e0)*HOLD;
		
		else if (HMIN>1.2E0 || (HMIN >= (double)(0.5) && HMIN<(double)0.9e0)) 
		
			H=HMIN*HOLD;
		
		else if (HMIN<(double)(0.5e0) ) 
		
			H=(double)(0.5e0)*HOLD;
		
		if (ERRORS>0) 
			ERRORS=0;
	}/*end IND=1*/


L27:	if (X<XEND){
		TEMP=XEND-X;
		
		if (TEMP<((double)(2.0e0*H))){
			if (TEMP>H) 
				H=TEMP*(double)(0.5e0);
				else 
				H=TEMP;
		}
				goto L10;
	
	}
g302: fprintf(f1,"TOL=%10.5e\n",TOL);
      fprintf(f1,"MAX ERROR=%10.5e\n",MAXE);
	  fprintf(f1,"XMAXE=%10.5e\n",XMAXE);
      fprintf(f1,"SUCCESS STEP=%d\n",SUCSTP);
      fprintf(f1,"FAIL STEP=%d\n",FA);
      fprintf(f1,"TOTAL STEPS=%d\n",TOTALSTP);
	  printf("TOL=%10.5e\n",TOL);
      printf("MAX ERROR=%10.5e\n",MAXE);
	  printf("XMAXE=%10.5e\n",XMAXE);
      printf("SUCCESS STEP=%d\n",SUCSTP);
      printf("FAIL STEP=%d\n",FA);
      printf("TOTAL STEPS=%d\n",TOTALSTP);
	fclose(f1);

	

}/*end INTEGRATE2*/


/*******************************************************************/
/*                        END OF SUBROUTINE INTEGRATE2             */
/*******************************************************************/
 
	/* BELOW ARE ALL THE SUBROUTINE CALL BY INTEGRATE2 */

/*******************************************************************/
/*                    EVALUATION OF FUNCTION                       */
/*******************************************************************/

void FCN(I,X,Y,F)
int I;
double X,Y[9][9],F[9];
   
 { 
	if(I==1)
	 F[1]=Y[1][0]*(4*pow(X,3)-Y[1][0])/(pow(X,4)-1);

	return;
	
  }

/*******************************************************************/
/*                   SUBROUTINE EXACT SOLUTIONS                    */
/*******************************************************************/

void SOLN()
 {
  
	int I;
	double SOL[9][9];

	SOL[1][0]=fabs(((1+X+pow(X,2)+pow(X,3))-Y[1][0])/(1+Y[1][0]));
	printf("SOL[%d][%d]=%15.10e\n",1,0,SOL[1][0]);			
    for(I=1;I<=N;I++)
	{
	  if(ERS[I]<SOL[I][0])
	  {
    	ERS[I]=SOL[I][0];
    	   if(MAXE<ERS[I])
		   {
	    	MAXE=ERS[I];
	    	XMAXE=X;
	    	HH=H;
		   }/*end inner if*/
	  }/*end outer if*/
	}/*end for loop*/
	return;
  } /* end of subroutine solutions*/

/*******************************************************************/
/*                           SUBROUTINE EVALUATE                   */
/*******************************************************************/

void EVALUATE(L,T,Y,F,R)
double Y[9][9],F[9];
int L,T,R; 
 {
    int SS4; 
    for (SS4=L;SS4<=T;SS4++)
	if ((R==0 && J[SS4]==0) || (R==1 && J[SS4]>0) || R==2)
	FCN(SS4,X,Y,F); /* call subroutine function*/
	return;
 }

/***************************************************************************/
/*                          SUBROUTINE STEPSIZE                            */
/***************************************************************************/

void STSIZE(J,D,Y,H)

int J[9],D[9];
double Y[9][9];
double *H;

{
	int I,MINS;
	double TEMP1;

      
    MINS=MAXD+1;
	MAX=(double)0.e0;
    
	for(I=1;I<=N;I++)
	{
	   if (MINS>S[I]+J[I]) 
		   MINS=S[I]+J[I];
		   TEMP1=(double)fabs(Y[I][D[I]]);
		
		if (MAX<TEMP1) 
			MAX=TEMP1;
	}/*end for loop*/

	if (MAX<1) 
		MAX=1;

	*H=(pow((TOL/MAX),(1/(double)(MINS+1))))/
           ((double)(pow(4.0,(double)(MINS))));
		
	if (*H>(XEND-X)/2)
		*H=*H/(double)(10);
return;
}

/***********************************************************************/
/*                          SUBROUTINE DERV                            */
/***********************************************************************/

void DERV(F1,Y,MJ,K)
int MJ[9],K[9];
double F1[9][14],Y[9][9];

{
	int I,FACT;
	
	for(I=1;I<=N;I++)
	{
		FACT=1;
		
		for(T=0;T<=K[I]-1;T++){

			F1[I][T]=Y[I][MJ[I]+1]/FACT;
			FACT=FACT*(T+1);
				}/*end T*/
		}/*end I*/
return;
	
}


the output should in the form of exponential. what can i do to make the output correct?

Is This A Good Question/Topic? 0
  • +

Replies To: the output of this program is SOL[1][0]=-1.#IND000000e+000.

#2 janotte  Icon User is offline

  • code > sword
  • member icon

Reputation: 991
  • View blog
  • Posts: 5,141
  • Joined: 28-September 06

Re: the output of this program is SOL[1][0]=-1.#IND000000e+000.

Posted 04 April 2011 - 03:56 AM

First step is to deal with these compiler warnings
dic.c:28: warning: return type defaults to ‘int’
dic.c: In function ‘main’:
dic.c:50: warning: control reaches end of non-void function
dic.c: In function ‘INTEGRATE2’:
dic.c:489: warning: label ‘L22’ defined but not used
dic.c:311: warning: label ‘L08’ defined but not used
dic.c:250: warning: label ‘L14’ defined but not used
dic.c:163: warning: label ‘L7’ defined but not used
dic.c:131: warning: label ‘L43’ defined but not used
dic.c:128: warning: label ‘L45’ defined but not used
dic.c:81: warning: unused variable ‘DEL’
dic.c:81: warning: unused variable ‘N12’
dic.c:81: warning: unused variable ‘N11’
dic.c:81: warning: unused variable ‘N10’
dic.c:79: warning: unused variable ‘DUMMY’
dic.c:77: warning: unused variable ‘EVALE2’
dic.c:77: warning: unused variable ‘FACT’
dic.c:76: warning: unused variable ‘IP1’
dic.c:74: warning: unused variable ‘ABSE’
dic.c:74: warning: unused variable ‘KH’
dic.c:74: warning: unused variable ‘MAXEG’
dic.c:72: warning: unused variable ‘TEMP4’
dic.c:72: warning: unused variable ‘TEMP3’
dic.c:72: warning: unused variable ‘E2’
dic.c:72: warning: unused variable ‘E1’
dic.c:71: warning: unused variable ‘KT’





Always enable the display of warnings in your compiler.
Error reporting is not enough to write good code.
If you don't know how to enable warnings in your tool then search the help files or the Google.
You need to learn how to display warnings as a matter of urgency. Don't write another line of code until you have added this vital skill to your toolbox.
Was This Post Helpful? 1
  • +
  • -

#3 nzm  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 3
  • Joined: 03-April 11

Re: the output of this program is SOL[1][0]=-1.#IND000000e+000.

Posted 07 April 2011 - 03:44 AM

For your information, this source code does not belong to me. So, it is kind of complicated and quite long. I'm just modified the code to solve some kind of higher order differential equations. This source code is also c programming code. I just want to know, is it the subroutine of evaluation of function is totally correct in line 582 to 594? The output of the previous code before i modified was not in the form SOL[1][0]=-1.#IND000000e+000. The output of the previous code is SOL[1][0]=5.3124999672e-009.

I have corrected the source code regarding the warnings.
The warning stated that FACT is unused variable but I used FACT in the SUBROUTINE DERV.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <usclkc.h>
#define true 1
#define false 0

/* GLOBAL VARIABLES DECLARATION*/

int N,TEST,MAXD,TESTS,STIFF;
double X,XEND,TOL,H;
double Y[9][9];
int D[9],J[9];

/*new add for global declaration*/

double HMIN,REL,AB,ED[9],RD[9],P[9][9],PSI[21],F[9],MAX,DELTA,YP[9],MAXE,XMAXE;
double F1[9][14],G[15][17],ER[9],JACBN[9][9],B[9],ERS[9],HH;
int FIRSTSTEP,STORE[9],FR,KM1,RR[9],K[9],S[9],NEWJAC,JUMP,CNT[9];
int ST[9],MJ[9],NS,MKP1,T,TT,JMAX,PS[9],R1; 

/**********************************************************************/
/*                         STARTS MAIN PROGRAM                        */
/**********************************************************************/

main()

{

	void INTEGRATE2();
	TOL=1.0E-7; /* TOL=tolerance limit */

	N=1;/* number of equations*/
	TEST=1;
	TOL=TOL*0.1;/* tolerance limit*/
	Y[1][0]=(double)15.e0;/* initial condition for test problem*/
	
	D[1]=1;/* order of the test problem*/
	

	X=(double)2.e0; /* initial value */
	XEND=(double)(25.e0);/* end point*/
	
	H=(double)1.e-4;
	J[1]=0;
	
    MAXD=1;

	INTEGRATE2(); /* Enter subroutine Integrate2*/
}

/*********************************************************************/
/*                     END OF MAIN PROGRAM                           */
/*********************************************************************/

/*********************************************************************/
/*                  STARTS SUBROUTINE INTEGRATE2                     */
/*********************************************************************/

void INTEGRATE2()

{

/* FUNCTIONS DECLARATION*/

void FCN(),SOLN(),STSIZE(),EVALUATE(),DERV();


/* VARIABLES DECLARATION*/

double TEMP,TEMP1,TEMP2,HOLD,XOLD,ERRMAX,EPS,
       RAT,SIGMA,SUM;

double E[16],DC[15][17],F2[9][4];
  
int    I,L,MAXJ,IM1,KP1,IND,KOLD,ERRORS,STEPS,
       K2,SS,GMAX,JP1,JP2,FACT;
 
int    SUCSTP,FA,TOTALSTP,BOOL[9];

double N13;

 FILE*f1;
	f1=fopen("C:PROFPRGRM.txt","w");

	if (TEST==2) REL=(double)0.e0; 
	else 
	    REL=(double)1.e0;

	if (TEST==3) AB=(double)0.e0; 
	else 
	    AB=(double)1.e0;
	
	if (TOL>1.0e-9) 
			FR=3; 
		else 
            FR=4;

        MAXE=(double)0.e0;
	    XMAXE=(double)0.e0;

	    HMIN = HOLD=(double)0.e0;
	    NS=0;

	 // TESTS=0; 
	    JMAX=0; 
	    EPS=1.0e-48; 
	    JUMP=false;
	    
		for (I=1;I<=N;I++)
		{
			BOOL[I]=false; 
	    	ERS[I]=(double)0.e0;
		}

	    for (I=1;I<=4;I++)
		    ER[I]=0.e0;
			ED[I]=0.e0;
		    
	    FIRSTSTEP=true;
	    STEPS=0;
		SUCSTP=0;
		FA=0;
		TOTALSTP=0;
	   	  
	    EVALUATE(1,N,Y,F,2);
	
/*L45:*/ 		for (I=1;I<=N;I++) 
			Y[I][D[I]]=F[I];

/*L43:*/ 		for (I=1;I<=N;I++)
			{
    			MJ[I]=D[I]-J[I];
     		    

		    if (MJ[I]==0) 
				S[I]=0;
			else if (MJ[I]==1) 
				S[I]=1;
			else if (MJ[I] <= 4) 
				S[I]=2; 
			else 
				S[I]=3;

		        RR[I]=J[I];

			}/*end I*/

	       if(H==0)STSIZE(J,D,Y,&H);
	       fprintf(f1,"Initial Stepsize=%10.5e\n\n",(double)H);
  
L4:     MKP1=0;
		for (I=1;I<=N;I++)
		{
			K[I]=J[I]+1;

			if (K[I]>MKP1) MKP1=K[I];
	    }/*end I*/

	    MKP1=MKP1+1;
    	DERV(F1,Y,MJ,K);
	    
/*L7:*/    	for(I=1;I<=MKP1;I++)
	
		PSI[I]=(double)0.e0;
		
/*********************************************************************/
/*                      COMPUTE CHANGE COEFFICIENTS                  */
/*********************************************************************/
L10:    if (H != HOLD) NS=0;
		
		if (NS >= MKP1) goto L11;
	
		if (NS<MKP1) NS=NS+1;

		if(MKP1>NS)
			for (I=MKP1;I>=NS+1;I--)
			{
			    PSI[I]=PSI[I-1]+H;
			}
				PSI[NS]=H*(double)(NS);
						
				for (I=NS;I<=MKP1;I++)
				{
		
					G[I][0]=PSI[I]*G[I-1][0];
				}	
				for (K2=1;K2<=N;K2++)
			
					if (MJ[K2] != 0)
					{
						GMAX=0;
						for (I=1;I<=N;I++)
						{
					    	L=K[I]+MJ[I];
						
							if (GMAX<L) 
								GMAX=L;
						}/*end I*/

					if (NS==1)
						{
							for (T=1;T<=GMAX+1;T++)
							{
								G[1][T-1]=G[0][T-1]*H/(double)(T);
								G[0][T]=G[1][T-1];
							}
						}/*end NS==1*/

					if ((MKP1>KOLD) && (NS>1))
						{
							G[1][GMAX]=H*G[1][GMAX-1]/(double)(GMAX+1);
							G[0][GMAX+1]=G[1][GMAX];

							for (I=2;I<=NS-1;I++)
							{
								T=GMAX+1-I;
								IM1=I-1;
								G[I][T]=PSI[I]*G[IM1][T]-((double)(T)*G[IM1][T+1]);
							}/*end I*/
						}/*end MKP1>KOLD*/

					if (NS==1) 
							L=2; 
						else 
							L=NS;

						for (I=L;I<=MKP1;I++)
						{
							IM1=I-1; 
							MAXJ= GMAX+1-I;
						    for (T=1;T<=MAXJ;T++)
							{
								G[I][T]=PSI[I]*G[IM1][T]-T*G[IM1][T+1];
							}
						}/*end I*/
						// goto L07;
					}/*end TT<N*/


		

L11:			XOLD=X;
				X=X+H;
				fprintf(f1,"X=%15.10e,\tH=%15.10e\n\n",X,H);
				printf("X=%15.10e,\tH=%15.10e\n\n",X,H);
				HOLD=H;
				KOLD=MKP1;

/*L14:*/  			HMIN=(double)3.e0;
				TOTALSTP=TOTALSTP+1;

	  
/************************************************************************/
/*                                  PREDICTION                          */
/************************************************************************/

for (I=1;I<=N;I++)
			 {
				 fprintf(f1,"K[%d]=%d\n\n",I,K[I]);
				 printf("K[%d]=%d\n\n",I,K[I]);
				for (K2=0;K2<=MJ[I];K2++)
				{
				  SIGMA=0.0;
				  if (K2==0) SUM=(double)0.e0;
				  else 
				  SUM=Y[I][MJ[I]-1];
				   for (T=0;T<=K[I]-1;T++)
						SIGMA=SIGMA + G[T][K2]*F1[I][T];
						  	if (K2<2)  goto L49;
												
							for (L=(K2-1);L>=1;L--)
							
								SUM=SUM*(H/L)+Y[I][MJ[I]-K2+L-1];

L49:	P[I][MJ[I]-K2]=SUM + SIGMA;/*predicted value*/
			fprintf(f1,"P[%d][%d]=%15.10e\n",I,MJ[I]-K2,P[I][MJ[I]-K2]);
				}/* end K2*/
							
						
					}/* end I*/


EVALUATE(1,N,P,F,0);
	for (I=1;I<=N;I++)
	if (J[I]==0)
	{
		ED[I]=F[I]-P[I][D[I]];
		for (L=1;L<=D[I];L++)
			P[I][MJ[I]-L]=P[I][MJ[I]-L]+ED[I]*G[K[I]][L]/G[K[I]][0];
	}
	    fprintf(f1,"P[%d][%d]=%15.10e\n",I,MJ[I]-L,P[I][MJ[I]-L]);						
		if(TT>=1) EVALUATE(1,N,P,F,1); /*try yang ni */


/***********************************************************************/
/*                                 ERROR COMPUTATION                   */
/***********************************************************************/


	ERRMAX=(double)0.e0;
	
	for(I=1;I<=N;I++){
		KM1=K[I]-1;

		TEMP=AB+REL*fabs(P[I][MJ[I]-S[I]]);
	if (J[I]==0)
	{
	E[I]=fabs((S[I])*G[KM1][S[I]+1]*ED[I]/G[K[I]][0])/TEMP;
	}
/*L08:*/	if (ERRMAX<E[I]){
			ERRMAX=E[I];
		
	}/*end I*/

	if (ERRMAX<TOL) 
	{
		IND=1;
	}
	else{
		  FA++;
		  IND=2;
		  goto L15;
	}

SUCSTP=SUCSTP+1; 
}
	
/***********************************************************************/
/*                                  CORRECTIONS                        */
/***********************************************************************/

	 EVALUATE(1,N,P,F,0);
/***************************************************************/
/*                    DIVIDED DIFFERENCES                      */
/***************************************************************/

L15:	if (IND==1){
		for (I=1;I<=N;I++){
	
			Y[I][D[I]]=F[I];

			for (L=1;L<=D[I];L++) 
				Y[I][D[I]-L]=P[I][D[I]-L];
			
			if (J[I]==0) 
				TEMP=F[I];
			else 
				TEMP=P[I][MJ[I]];

			if (FIRSTSTEP) 
				T=K[I]-1;
			else 
				T=K[I];

			for (L=0;L<=T;L++){
				TEMP1=(TEMP-F1[I][L])/PSI[L+1];
				F1[I][L]=TEMP; 
				TEMP=TEMP1;
			}/*end for statement L*/
			
			F1[I][T+1]=TEMP;
		}/*end I*/
		SOLN();
	}/*end IND==1*/
	else{	
		for (I=1;I<=N;I++){
			F2[I][3]=ED[I]/(G[K[I]-1][0]*PSI[K[I]]);

			if (K[I]>J[I]+1)
				F2[I][2]=F1[I][K[I]-1]+PSI[K[I]]*F2[I][3];

			if (K[I]>J[I]+2)
				F2[I][1]=F1[I][K[I]-2]+PSI[K[I]-1]*F2[I][2]; 
			
		}/*end IND==2*/
	}/*end I*/

if (FIRSTSTEP){
	        if (IND==1)
			{
		    	FIRSTSTEP=false;
			    for (I=1;I<=N;I++)
					ERRORS=0;
	           goto L27;
			}
				else 
		{
			if (ERRORS==0) 
				HMIN=(double)1.e0/(double)(1.e1);
			else if (ERRORS==1)
				HMIN=(double)1.e0/(double)(1.e3);
		}
				}/*end FIRSTSTEP*/

/**************************************************************************/
/*                              ORDER AND STEPSIZE                        */
/**************************************************************************/

MKP1=0;
		for (I=1;I<=N;I++)
	{	
		KP1=K[I]+1;
		JP1=J[I]+1;
		KM1=K[I]-1;
		JP2=JP1+1;

       if (K[I]>JP1+1)
	   	      SS=1;
	   else if (K[I]>JP1){
		   SS=2;}
	   else {
		   SS=3;}
   		for(L=SS;L<=3;L++)
		{
		T=KM1+L-2;
		if (IND==1)
			TEMP=F1[I][T];
			else 
			  TEMP=F2[I][L];
			if (J[I]==0)
				ER[L]=(double)fabs(S[I]*G[T-1][S[I]+1]*TEMP);
			
		}/*end L*/
		if ((K[I]==JP2) && (ER[2]<(0.5*ER[3])))
			 K[I]=KM1;
			
		if (K[I]>JP2)
			{
				if (ER[1]>ER[2])
					TEMP1=ER[1];
				else 
					TEMP1=ER[2];
				if (TEMP1<ER[3]) 
			    	K[I]=KM1;
			
			}
		if(J[I]==0)
			N13=K[I]<12;
			else 
			N13=K[I]<6;
		if ((! BOOL [I]) && (K[I] != KM1) && (IND==1) && (N13))
		{
			RAT=(double)(2.e0);
			K[I]=KP1;
			goto L12;
		} /*end BOOL[I]*/
		else 
			BOOL[I]=true;
		if (IND==1){
			if((NS >= KP1) && (K[I] != KM1) && (J[I]==0 && K[I]<12) /*|| (J[I]>0 && K[I]<6)*/) 
{
  if (J[I]==0) 
	ER[4]=(double)fabs(S[I]*G[K[I]][S[I]+1]*F1[I][KP1]);
	else 
    ER[4]=(double)fabs((DC[K[I]][J[I]]/DC[KP1][J[I]]*G[KP1][S[I]]-G[K[I]][S[I]])*PSI[KP1]*F1[I][KP1]);
		if (K[I]==JP1)
		{
		 if (ER[4]<(double)(0.5)*ER[3]) 
			K[I]=KP1;
		}
		else
		{
		  	if ((/*(J[I]>0 && ER[4]<ER[3]) || */ (J[I]==0 && ER[4]<ER[3])) && ER[3]<ER[2])
			K[I]=KP1;
		    else if (ER[2]<ER[3] && ER[2]<ER[4])
		        	K[I]=KM1;
		}/*end K[I]>1*/
	}/*end IND=1*/

}/*end IND=1*/

	TEMP2=AB+REL*(double)fabs(P[I][MJ[I]-S[I]]);
	T=K[I]-KM1+2;

	if (ER[T]<EPS) 
		RAT=(double)(4.e0);
	else 
		RAT=(double)pow(TOL*TEMP2/ER[T],(double)1.e0/(double)(K[I]+S[I]));
	

L12:	if (HMIN>RAT) HMIN=RAT;
		if (MKP1<K[I]) MKP1=K[I];

	}/*end I*/
	
	MKP1=MKP1+1;

/*L22:*/	if (IND==2){
		ERRORS=ERRORS+1;
	    
		if (! FIRSTSTEP){
			if (HMIN <=(double)0.5e0) 
				HMIN=0.4;
			else if (HMIN>(double)0.5e0 && HMIN <=(double)0.8e0) 
				HMIN=(double)0.5e0;
				else 
				HMIN=(double)0.7e0;
				}

		H=HMIN*HOLD;
		X=XOLD;

		if (NS==0) 
			T=2;
		else 
			T=NS+1;
	
		for(T=NS+1;T<=KOLD;T++) 
			PSI[T-1]=PSI[T]-HOLD;
	 

		if (ERRORS==12){
			printf("TOO MANY ERRORS %20.5e %d %d\n",X,0,5);
			goto g302;
		}
		else if (ERRORS>3){
			H=H/(double)(10);
			goto L4;
		}
	}/*end IND==2*/
	else{
		HMIN=HMIN*(double)0.8e0;
		
		if (HMIN >=(double)(2)) 
		
			H=(double)(2.e0)*HOLD;
		
		else if (HMIN>1.2E0 || (HMIN >= (double)(0.5) && HMIN<(double)0.9e0)) 
		
			H=HMIN*HOLD;
		
		else if (HMIN<(double)(0.5e0) ) 
		
			H=(double)(0.5e0)*HOLD;
		
		if (ERRORS>0) 
			ERRORS=0;
	}/*end IND=1*/


L27:	if (X<XEND){
		TEMP=XEND-X;
		
		if (TEMP<((double)(2.0e0*H))){
			if (TEMP>H) 
				H=TEMP*(double)(0.5e0);
				else 
				H=TEMP;
		}
				goto L10;
	
	}
g302: fprintf(f1,"TOL=%10.5e\n",TOL);
      fprintf(f1,"MAX ERROR=%10.5e\n",MAXE);
	  fprintf(f1,"XMAXE=%10.5e\n",XMAXE);
      fprintf(f1,"SUCCESS STEP=%d\n",SUCSTP);
      fprintf(f1,"FAIL STEP=%d\n",FA);
      fprintf(f1,"TOTAL STEPS=%d\n",TOTALSTP);
	  printf("TOL=%10.5e\n",TOL);
      printf("MAX ERROR=%10.5e\n",MAXE);
	  printf("XMAXE=%10.5e\n",XMAXE);
      printf("SUCCESS STEP=%d\n",SUCSTP);
      printf("FAIL STEP=%d\n",FA);
      printf("TOTAL STEPS=%d\n",TOTALSTP);
	fclose(f1);

	

}/*end INTEGRATE2*/


/*******************************************************************/
/*                        END OF SUBROUTINE INTEGRATE2             */
/*******************************************************************/
 
	/* BELOW ARE ALL THE SUBROUTINE CALL BY INTEGRATE2 */

/*******************************************************************/
/*                    EVALUATION OF FUNCTION                       */
/*******************************************************************/

void FCN(I,X,Y,F)
int I;
double X,Y[9][9],F[9];
   
 { 
	if(I==1)
	 F[1]=Y[1][0]*(4*pow(X,3)-Y[1][0])/(pow(X,4)-1);

	return;
	
  }

/*******************************************************************/
/*                   SUBROUTINE EXACT SOLUTIONS                    */
/*******************************************************************/

void SOLN()
 {
  
	int I;
	double SOL[9][9];

	SOL[1][0]=fabs(((1+X+pow(X,2)+pow(X,3))-Y[1][0])/(1+Y[1][0]));
	printf("SOL[%d][%d]=%15.10e\n",1,0,SOL[1][0]);			
    for(I=1;I<=N;I++)
	{
	  if(ERS[I]<SOL[I][0])
	  {
    	ERS[I]=SOL[I][0];
    	   if(MAXE<ERS[I])
		   {
	    	MAXE=ERS[I];
	    	XMAXE=X;
	    	HH=H;
		   }/*end inner if*/
	  }/*end outer if*/
	}/*end for loop*/
	return;
  } /* end of subroutine solutions*/

/*******************************************************************/
/*                           SUBROUTINE EVALUATE                   */
/*******************************************************************/

void EVALUATE(L,T,Y,F,R)
double Y[9][9],F[9];
int L,T,R; 
 {
    int SS4; 
    for (SS4=L;SS4<=T;SS4++)
	if ((R==0 && J[SS4]==0) || (R==1 && J[SS4]>0) || R==2)
	FCN(SS4,X,Y,F); /* call subroutine function*/
	return;
 }

/***************************************************************************/
/*                          SUBROUTINE STEPSIZE                            */
/***************************************************************************/

void STSIZE(J,D,Y,H)

int J[9],D[9];
double Y[9][9];
double *H;

{
	int I,MINS;
	double TEMP1;

      
    MINS=MAXD+1;
	MAX=(double)0.e0;
    
	for(I=1;I<=N;I++)
	{
	   if (MINS>S[I]+J[I]) 
		   MINS=S[I]+J[I];
		   TEMP1=(double)fabs(Y[I][D[I]]);
		
		if (MAX<TEMP1) 
			MAX=TEMP1;
	}/*end for loop*/

	if (MAX<1) 
		MAX=1;

	*H=(pow((TOL/MAX),(1/(double)(MINS+1))))/
           ((double)(pow(4.0,(double)(MINS))));
		
	if (*H>(XEND-X)/2)
		*H=*H/(double)(10);
return;
}

/***********************************************************************/
/*                          SUBROUTINE DERV                            */
/***********************************************************************/

void DERV(F1,Y,MJ,K)
int MJ[9],K[9];
double F1[9][14],Y[9][9];

{
	int I,FACT;
	
	for(I=1;I<=N;I++)
	{
		FACT=1;
		
		for(T=0;T<=K[I]-1;T++){

			F1[I][T]=Y[I][MJ[I]+1]/FACT;
			FACT=FACT*(T+1);
				}/*end T*/
		}/*end I*/
return;
	
}

This post has been edited by nzm: 07 April 2011 - 03:50 AM

Was This Post Helpful? 0
  • +
  • -

#4 JackOfAllTrades  Icon User is offline

  • Saucy!
  • member icon

Reputation: 6246
  • View blog
  • Posts: 24,014
  • Joined: 23-August 08

Re: the output of this program is SOL[1][0]=-1.#IND000000e+000.

Posted 07 April 2011 - 04:03 AM

No. You have a GLOBAL FACT variable which you never use. You use a LOCAL FACT variable in the DERV function.

Read about functions in the tutorials linked in my signature.

Stylistically this code is horrible for the following reasons:

1. too many global variables
2. using UPPERCASE for all variable names
3. poor variable naming
4. use of an implicit main -- should be int main(void)
5. use of goto

People just aren't going to want to help with such horrible code. I hurts our heads.
Was This Post Helpful? 0
  • +
  • -

#5 nzm  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 3
  • Joined: 03-April 11

Re: the output of this program is SOL[1][0]=-1.#IND000000e+000.

Posted 08 April 2011 - 01:33 AM

But, what can I do if the output in that form?
What is actually the reason that the output could become like that?
I have never seen the output in that form before.
I know it is really hard for anyone to revise and help me regarding the horrible code, but I really need help to find out what is the reason of the output could become like that.
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1