/*********************************************************************** * Solve a boundary value problem for a first order DE system of the * * form: * * * * Y1' = F1(X,Y1,Y2,...,YN) * * Y2' = F2(X,Y1,Y2,...,YN) * * ... * * YN' = FN(X,Y1,Y2,...,YN) * * * * via the shooting method by determining an approximation for the * * initial value Y(A). * * -------------------------------------------------------------------- * * SAMPLE RUNS; * * * * TEST EXAMPLE (WITH METHOD #1): * * Method: Runge-Kutta embedding formula of 4/5th order. * * ============= * * ANALYZED SET OF DIFFERENTIAL EQUATIONS: * * --------------------------------------- * * Y'( 1) = Y(2) * * Y'( 2) = -Y(1)**3 * * WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * * AT POINT A = 0.000E+00 : 0.000E+00 * * AT POINT B = 1.000E+00 : 0.000E+00 * * * * AT LOCATION A = 0.000E+00 THE FOLLOWING VALUES ARE PROVIDED: * * * * Y( 1) = 0.000E+00 * * Y( 2) = 1.200E+01 * * * * REQUIRED PARAMETER: * * ------------------: * * -PRECISION FOR THE IVP = 1.0000000000E-07 * * -PRECISION FOR THE BVP = 1.0000000000E-08 * * -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * * -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * * -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * * * * SOLUTION OF THE PROBLEM: * * ------------------------ * * THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * * IN POINT A = 0.000E+00 AS FOLLOWS: * * * * Y( 1) = 0.0000000000E+000 * * Y( 2) = 9.7229810350E+000 * * * * THIS REQUIRES 3 NEWTON-ITERATIONS! * * * * DETERMINATION COMPLETED AS PLANNED! * * * * * * TEST EXAMPLE (WITH METHOD #2): * * Method: Predictor-corrector method of order 4 by Adams-Bashforth- * * Moulton. * * ============= * * ANALYZED SET OF DIFFERENTIAL EQUATIONS: * * --------------------------------------- * * Y'( 1) = Y(2) * * Y'( 2) = -Y(1)**3 * * WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * * AT POINT A = 0.000E+00 : 0.000E+00 * * AT POINT B = 1.000E+00 : 0.000E+00 * * * * AT LOCATION A = 0.000E+00 THE FOLLOWING VALUES ARE PROVIDED: * * * * Y( 1) = 0.000E+00 * * Y( 2) = 1.200E+01 * * * * REQUIRED PARAMETER: * * ------------------: * * -PRECISION FOR THE IVP = 1.0000000000E-07 * * -PRECISION FOR THE BVP = 1.0000000000E-08 * * -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * * -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * * -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * * * * SOLUTION OF THE PROBLEM: * * ------------------------ * * THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * * IN POINT A = 0.000E+00 AS FOLLOWS: * * * * Y( 1) = 0.0000000000E+000 * * Y( 2) = 9.7229809582E+000 * * * * THIS REQUIRES 4 NEWTON-ITERATIONS! * * * * DETERMINATION COMPLETED AS PLANNED! * * * * * * TEST EXAMPLE (WITH METHOD #3): * * Method: Runge-Kutta embedding formula of 7/8th order. * * ============= * * ANALYZED SET OF DIFFERENTIAL EQUATIONS: * * --------------------------------------- * * Y'( 1) = Y(2) * * Y'( 2) = -Y(1)**3 * * WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * * AT POINT A = 0.000E+00 : 0.000E+00 * * AT POINT B = 1.000E+00 : 0.000E+00 * * * * AT LOCATION A = 0.000E+00 THE FOLLOWING VALUES ARE PROVIDED: * * * * Y( 1) = 0.000E+00 * * Y( 2) = 1.200E+01 * * * * REQUIRED PARAMETER: * * ------------------: * * -PRECISION FOR THE IVP = 1.0000000000E-07 * * -PRECISION FOR THE BVP = 1.0000000000E-08 * * -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * * -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * * -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * * * * SOLUTION OF THE PROBLEM: * * ------------------------ * * THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * * IN POINT A = 0.000E+00 AS FOLLOWS: * * * * Y( 1) = 0.0000000000E+000 * * Y( 2) = 9.7229804225E+000 * * * * THIS REQUIRES 5 NEWTON-ITERATIONS! * * * * DETERMINATION COMPLETED AS PLANNED! * * * * * * TEST EXAMPLE (WITH METHOD #4): * * Method: Extrapolation method of Bulirsch-Stoer. * * ============= * * ANALYZED SET OF DIFFERENTIAL EQUATIONS: * * --------------------------------------- * * Y'( 1) = Y(2) * * Y'( 2) = -Y(1)**3 * * WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * * AT POINT A = 0.000E+00 : 0.000E+00 * * AT POINT B = 1.000E+00 : 0.000E+00 * * * * AT LOCATION A = 0.000E+00 THE FOLLOWING VALUES ARE PROVIDED: * * * * Y( 1) = 0.000E+00 * * Y( 2) = 1.200E+01 * * * * REQUIRED PARAMETER: * * ------------------: * * -PRECISION FOR THE IVP = 1.0000000000E-07 * * -PRECISION FOR THE BVP = 1.0000000000E-08 * * -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * * -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * * -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * * * * SOLUTION OF THE PROBLEM: * * ------------------------ * * THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * * IN POINT A = 0.000E+00 AS FOLLOWS: * * * * Y( 1) = 0.0000000000E+000 * * Y( 2) = 9.7229810093E+000 * * * * THIS REQUIRES 3 NEWTON-ITERATIONS! * * * * DETERMINATION COMPLETED AS PLANNED! * * * * * * TEST EXAMPLE (WITH METHOD #5): * * Method: Implicit Runge-Kutta-Gauss method. * * ============= * * ANALYZED SET OF DIFFERENTIAL EQUATIONS: * * --------------------------------------- * * Y'( 1) = Y(2) * * Y'( 2) = -Y(1)**3 * * WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * * AT POINT A = 0.000E+00 : 0.000E+00 * * AT POINT B = 1.000E+00 : 0.000E+00 * * * * AT LOCATION A = 0.000E+00 THE FOLLOWING VALUES ARE PROVIDED: * * * * Y( 1) = 0.000E+00 * * Y( 2) = 1.200E+01 * * * * REQUIRED PARAMETER: * * ------------------: * * -PRECISION FOR THE IVP = 1.0000000000E-07 * * -PRECISION FOR THE BVP = 1.0000000000E-08 * * -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * * -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * * -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * * * * SOLUTION OF THE PROBLEM: * * ------------------------ * * THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * * IN POINT A = 0.000E+00 AS FOLLOWS: * * * * Y( 1) = 0.0000000000E+000 * * Y( 2) = 9.7229810243E+000 * * * * THIS REQUIRES 4 NEWTON-ITERATIONS! * * * * DETERMINATION COMPLETED AS PLANNED! * * * * -------------------------------------------------------------------- * * REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * * and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * * * * Visual C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ // Initialize data; test examples can be exchanged by adjusting // the functions dgl(), randbed() and the macro N. #include // for umleiten, fprintf, stderr, printf, // scanf, NULL, REAL, LZS, LZP, // fehler_melden #include // for vmalloc, vmcomplete, vmfree, // vminit, VEKTOR #include "rwp.h" // for rwp // table of method names static char *TMethod[]= { "Runge-Kutta embedding formula of 4/5th order.", "Predictor-corrector method of order 4 by Adams-Bashforth-Moulton.", "Runge-Kutta embedding formula of 7/8th order.", "Extrapolation method of Bulirsch-Stoer.", "Implicit Runge-Kutta-Gauss method.", }; void Pause() { char s[2]; printf("Any key + Enter to start..."); scanf("%s",s); } static void dgl(REAL x, REAL *y, REAL *f) { x = x; f[0] = y[1]; f[1] = -y[0] * y[0] * y[0]; } static void randbed(REAL *ya, REAL *yb, REAL *r) { r[0] = ya[0]; r[1] = yb[0]; } int main(void) { #define N 2 // Number of equations static REAL yanf0[N] = { ZERO, (REAL)12.0 }; static char ystri[N][80] = { "Y(2)", // alpha-numeric form "-Y(1)**3" }; // of the DE REAL a = ZERO; REAL b = ONE; REAL h = (REAL)0.01; REAL wa = ZERO; REAL wb = ZERO; long ifmax = 20000; int itmax = 1000; REAL epsawp = (REAL)1.0e-7; REAL epsrb = (REAL)1.0e-8; REAL yanf[N]; int fehler; int iter; int method; int j; #define print(text) printf("* %s\n", text); // output of the test examples for five varying methods Pause(); for (method = 1; method <= 5; method++) { copy_vector(yanf, yanf0, N); // yanf = yanf0 print(""); print(""); print(""); printf("* TEST EXAMPLE (WITH METHOD #%1d):\n",method); printf("* Method: %s\n", TMethod[method-1]); print("============="); print("ANALYZED SET OF DIFFERENTIAL EQUATIONS:"); print("---------------------------------------"); for (j = 0; j < N; j++) printf("* Y'(%2d) = %s\n", j + 1, ystri[j]); print("WITH BOUNDARY CONDITIONS : VALUE OF Y(1)"); printf("* AT POINT A = %10.3"LZP"E :%10.3"LZP"E\n" "* AT POINT B = %10.3"LZP"E :%10.3"LZP"E\n", a, wa, b, wb); print(""); printf("* AT LOCATION A = %10.3"LZP"E THE FOLLOWING " "VALUES ARE PROVIDED:\n", a); print(""); for (j = 0; j < N; j++) printf("* Y(%2d) = %10.3"LZP"E\n", j + 1, yanf[j]); print(""); print("REQUIRED PARAMETER:"); print("------------------:"); printf("* -PRECISION FOR THE IVP =%20.10"LZP"E\n" "* -PRECISION FOR THE BVP =%20.10"LZP"E\n" "* -STEP WIDTH FOR THE IVP=%20.10"LZP"E\n" "* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE " "IVP'S = %6ld\n" "* -MAX. NUMBER OF NEWTON-ITERATION STEPS " " = %6d\n", epsawp, epsrb, h, ifmax, itmax ); fehler = rwp(a, b, h, yanf, N, dgl, randbed, method, epsawp, epsrb, ifmax, itmax, &iter); // report result if (fehler != 0) { print(""); switch (fehler) { case 1: print("ERRORBOUND(S) TOO SMALL"); break; case 2: print("ERROR: B <= A"); break; case 3: print("ERROR: STEP SIZE H <= 0"); break; case 4: print("ERROR: NUMBER N OF DEQ'S NOT CORRECT"); break; case 5: print("ERROR: WRONG IVP-PROGRAM CHOICE"); break; case 6: print("IFMAX TO SMALL FOR SOLUTION OF THE " "IVP-PROGRAM"); break; case 7: print("AFTER ITMAX STEPS PRECISION WAS NOT REACHED!"); break; case 8: print("ERROR: JACOBI-MATRIX IS SINGULAR!"); break; case 9: print("LACK of MEMORY!"); break; default: print("UNKNOWN ERROR"); printf("* ERROR NUMBER = %2d\n", fehler); break; } } else { // no error? // Output of solution print(""); print("SOLUTION OF THE PROBLEM:"); print("------------------------"); print("THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED"); printf("* IN POINT A =%10.3"LZP"E AS FOLLOWS:\n", a); print(""); for (j = 0; j < N; j++) printf("* Y(%2d) = %20.10"LZP"E\n", j + 1, yanf[j]); print(""); printf("* THIS REQUIRES %5d NEWTON-ITERATIONS!\n", iter); print(""); print("DETERMINATION COMPLETED AS PLANNED!"); } Pause(); } return 0; } // -------------------------- END m_rwp.cpp ----------------------------