/******************************************************************* * Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the * * Euler-Romberg Method * * ---------------------------------------------------------------- * * REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en * * Basic et en Pascal By Claude Nowakowski, Edition du * * P.S.I., 1984" [BIBLI 04]. * * ---------------------------------------------------------------- * * SAMPLE RUN: * * (Solve Y' = X*X*Y with Y(0)=1 for X0=0 up to X1=1.0, exact * * solution is Y = Exp(X^3/3) ). * * * * Starting X...: 0 * * Ending X.....: 1 * * Value Y at X0: 1 * * Initial step : 0.1 * * Desired error: 1e-6 * * * * X Y Y Error Number of * * estimated exact subdivisions * * ---------------------------------------------------------- * * 0.1 1.000333 1.000333 0.00000000 4 * * 0.2 1.002670 1.002670 0.00000001 4 * * 0.3 1.009041 1.009041 0.00000006 4 * * 0.4 1.021562 1.021562 0.00000014 4 * * 0.5 1.042547 1.042547 0.00000027 4 * * 0.6 1.074654 1.074655 0.00000086 4 * * 0.7 1.121125 1.121126 0.00000107 4 * * 0.8 1.186094 1.186095 0.00000126 4 * * 0.9 1.275067 1.275069 0.00000133 4 * * 1.0 1.395611 1.395612 0.00000114 4 * * 1.1 1.558410 1.558411 0.00000047 4 * * ---------------------------------------------------------- * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *******************************************************************/ #include #include const DIM = 100; const LA=10; // maximum number of subdivisions double X[DIM], Y[DIM]; int NL[DIM]; int N,NC; double EF,ER,H,X0,X1,Y0,YEX; // Y' = F(X,Y) double F(double X, double Y) { return X*X*Y; } // Exact solution FX(X) double FX(double X) { return exp(X*X*X/3); } void Euler_Romberg(int NC, double H, double ER, double X0, double Y0, double *X, double *Y, int *NL) { /************************************************************** * Solve an Ordinary Differential Equation Y'=F(X,Y) using the * * Euler-Romberg Method. * * ----------------------------------------------------------- * * Inputs: * * NC: number of points to calculate * * H: abscissa integration step * * ER: desired precision * * X0: starting abscissa * * Y0: initial value of Y at X=X0 * * F: external user defined function Y' * * Outputs: * * X: Table of NC abscissas * * Y: Table of NC solution ordinates * * NL: Table of number of steps per abscissa * * * * Note: the ending abscissa X1 is given by X0 + NC*H. * **************************************************************/ int N; double ET,XC,YC; double T[20]; int J,K,L,LM,M,MM; // Initial conditions X[0]=X0; Y[0]=Y0; //main integration loop for (N=0; N<=NC; N++) { XC=X[N]; YC=Y[N]; T[1]=Y[N] + H*F(XC,YC); L=1; LM=2; do { XC=X[N]; YC=Y[N]; for (J=1; J<=LM; J++) { XC=XC+H/LM; YC=YC+H/LM*F(XC,YC); } T[L+1]=YC; M=1; K=L; MM=2; ET=1.0; if (K>1) do { T[K]=(MM*T[K+1]-T[K])/(MM-1); ET=fabs(T[K]-T[K-1]); M++; K--; MM *= 2; } while (ET>=ER && K>1); if (K==1) { L++; LM *= 2; } } while (L=ER); X[N+1]=X[N]+H; Y[N+1]=T[K]; } } void main() { // Initial data printf("\n"); printf(" Starting X...: "); scanf("%lf", &X0); printf(" Ending X.....: "); scanf("%lf", &X1); printf(" Value Y at X0: "); scanf("%lf", &Y0); printf(" Initial step : "); scanf("%lf", &H); printf(" Desired error: "); scanf("%lf", &ER); //number of calculations NC=(int) ((X1-X0)/H); //call integration routine Euler_Romberg(NC,H,ER,X0,Y0,X,Y,NL); //write header printf("\n"); printf(" X Y Y Error Number of \n"); printf(" estimated exact subdivisions \n"); printf("---------------------------------------------------------\n"); //results loop for (N=0; N<=NC; N++) { YEX=FX(X[N+1]); EF=fabs(YEX-Y[N+1]); printf("%6.2f%12.6f%12.6f%12.8f%8d\n", X[N+1],Y[N+1],YEX,EF,NL[N+1]); } printf("---------------------------------------------------------\n"); } // end of file eulromb.cpp