{*********************************************************************** * * * Test examples for methods to solve first order ordinary systems of * * differential equations. * * * * We supply the right hand sides, their alpha-numeric description and * * the exact analytic solution, if known. * * * * When running the main test program, the user can select among the * * examples below. * * * * TPW Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * -------------------------------------------------------------------- * * REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * * and Frank Uhlig, Springer-Verlag, 1996". * ***********************************************************************} UNIT t_dgls; INTERFACE Uses FGauss; Var bspnummer, n: Integer; dgltxt: Array[0..4] of String[70]; Procedure dgl(x:double; y:pVEC; f:pVEC); IMPLEMENTATION {Note: here, only examples #0 to #4 are implemented in Pascal} Procedure dgl(x:double; y:pVEC; f:pVEC); Begin Case bspnummer of 0: begin n:=2; f^[0] := y^[0] * y^[1] + COS(x) - 0.5 * SIN(2.0 * x); f^[1] := y^[0] * y^[0] + y^[1] * y^[1] - (ONE + SIN(x)); dgltxt[0]:=' y1'' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)'; dgltxt[1]:=' y2'' = y1 * y1 + y2 * y2 - (1 + sin(x))' end; 1: begin n:=1; f^[0]:=-Y^[0] + x/((1+x)*(1+x)); dgltxt[0]:=' y'' = -y + x/((1+x)*(1+x))' end; 2: begin {example for m_rwp} n:=2; f^[0] := y^[1]; f^[1] := -y^[0] * y^[0] * y^[0]; dgltxt[0]:=' y1'' = y2'; dgltxt[1]:=' y2'' = -y1^3' end; 3: begin n:=5; f^[0] := y^[1]; f^[1] := y^[2]; f^[2] := y^[3]; f^[3] := y^[4]; f^[4] := (45.0 * y^[2] * y^[3] * y^[4] - 40.0 * y^[3] * y^[3] * y^[3]) / (9.0 * y^[2] * y^[2]); dgltxt[0]:=' y1'' = y2'; dgltxt[1]:=' y2'' = y3'; dgltxt[2]:=' y3'' = y4'; dgltxt[3]:=' y4'' = y5'; dgltxt[4]:=' y5'' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)' end; 4: begin {y"=(2y'y'+yy)/y} n:=2; f^[0] := y^[1]; f^[1] := (2.0*y^[1]*y^[1]+y^[0]*y^[0])/y^[0]; dgltxt[0]:=' y1'' = y2'; dgltxt[1]:=' y2'' = (2y2*y2+y1*y1)/y1' end end End; END. {dgl} {Other examples given in C} (* ----------------------- DE system number 0 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl0(REAL x, REAL *y, REAL *f^) { f^[0] = y^[0] * y^[1] + COS(x) - HALf^ * SIN(TWO * x); f^[1] = y^[0] * y^[0] + y^[1] * y^[1] - (ONE + SIN(x)); } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt0(void) { return "y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)\n" "y2' = y1 * y1 + y2 * y2 - (1 + sin(x))\n"; } /* ----------------------- DE system number 1 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl1(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; f^[1] = (ONE + y^[1] * y^[1]) / y^[0]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt1(void) { return "y1' = y2\n" "y2' = (1 + y2 * y2) / y1\n"; } /* ----------------------- DE system number 2 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl2(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; f^[1] = -f^OUR * y^[0] + EXP(x); } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt2(void) { return "y1' = y2\n" "y2' = -4 * y1 + exp(x)\n"; } /* ----------------------- DE system number 3 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl3(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; f^[1] = y^[2]; f^[2] = y^[3]; f^[3] = y^[4]; f^[4] = ((REAL)45.0 * y^[2] * y^[3] * y^[4] - (REAL)40.0 * y^[3] * y^[3] * y^[3]) / (NINE * y^[2] * y^[2]); } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt3(void) { return "y1' = y2\n" "y2' = y3\n" "y3' = y4\n" "y4' = y5\n" "y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)\n"; } /* ----------------------- DE system number 4 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl4(REAL x, REAL *y, REAL *f^) { f^[0] = -y^[1]; f^[1] = y^[0]; f^[2] = y^[4] - y^[1] * y^[2]; f^[3] = y^[0] * y^[3] - y^[5]; f^[4] = -y^[2] - y^[1] * y^[4]; f^[5] = y^[3] + y^[0] * y^[5]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt4(void) { return "y1' = -y2\n" "y2' = y1\n" "y3' = y5 - y2 * y3\n" "y4' = y1 * y4 - y6\n" "y5' = -y3 - y2 * y5\n" "y6' = y4 + y1 * y6\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem y(0) = (1,0,0,1,e,0) * ***********************************************************************/ static void loesung4(REAL x, REAL *y) { y^[0] = COS(x); y^[1] = SIN(x); y^[2] = SIN(x) * EXP(COS(x)); y^[3] = COS(x) * EXP(SIN(x)); y^[4] = COS(x) * EXP(COS(x)); y^[5] = SIN(x) * EXP(SIN(x)); } /* ----------------------- DE system number 5 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl5(REAL x, REAL *y, REAL *f^) { f^[0] = (REAL)-500.5 * y^[0] + (REAL)499.5 * y^[1]; f^[1] = (REAL)499.5 * y^[0] - (REAL)500.5 * y^[1]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt5(void) { return "y1' = -500.5 * y1 + 499.5 * y2\n" "y2' = 499.5 * y1 - 500.5 * y2\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem y(0) = (4,2) * ***********************************************************************/ static void loesung5(REAL x, REAL *y) { y^[0] = THREE * EXP(-x) + EXP((REAL)-1000.0 * x); y^[1] = THREE * EXP(-x) - EXP((REAL)-1000.0 * x); } /* ----------------------- DE system number 6 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl6(REAL x, REAL *y, REAL *f^) { f^[0] = (REAL)15.0 * y^[1]; f^[1] = (REAL)-0.6 * y^[0]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt6(void) { return "y1' = 15.0 * y2\n" "y2' = -0.6 * y1\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem y(0) = (0,1) * ***********************************************************************/ static void loesung6(REAL x, REAL *y) { y^[0] = f^IVE * SIN(THREE * x); y^[1] = COS(THREE * x); } /* ----------------------- DE system number 7 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl7(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; f^[1] = y^[2]; f^[2] = y^[3]; f^[3] = f^OUR * y^[3] - SIX * y^[2] + f^OUR * y^[1] - y^[0]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt7(void) { return "y1' = y2\n" "y2' = y3\n" "y3' = y4\n" "y4' = 4.0 * y4 - 6.0 * y3 + 4.0 * y2 - y1\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem y(0) = (0,0,0,6) * ***********************************************************************/ static void loesung7(REAL x, REAL *y) { y^[0] = x * x * x * EXP(x); y^[1] = x * x * (x + THREE) * EXP(x); y^[2] = x * (x * x + SIX * x + SIX) * EXP(x); y^[3] = (x * x * x + NINE * x * x + (REAL)18.0 * x + SIX) * EXP(x); } /* ----------------------- DE system number 8 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl8(REAL x, REAL *y, REAL *f^) { f^[0] = (REAL)-52.0 * y^[0] + (REAL)50.0 * y^[1]; f^[1] = (REAL)50.0 * y^[0] - (REAL)52.0 * y^[1]; } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt8(void) { return "y1' = -52 * y1 + 50 * y2\n" "y2' = 50 * y1 - 52 * y2\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem * * y(0) = (1265588.55375228,-1265586.55375225) * ***********************************************************************/ static void loesung8(REAL x, REAL *y) { y^[0] = (REAL)1.000000015 * EXP( -TWO * x) + (REAL)1265587.553752265 * EXP((REAL)-102.0 * x); y^[1] = (REAL)1.000000015 * EXP( -TWO * x) - (REAL)1265587.553752265 * EXP((REAL)-102.0 * x); } /* ----------------------- DE system number 9 ----------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl9(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; #if^ def^ined(BC3) _f^preset(); #endif^ f^[1] = -y^[0] * COSH(x); } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt9(void) { return "y1' = y2\n" "y2' = -y1 * cosh(x)\n"; } /* ----------------------- DE system number 10 --------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl10(REAL x, REAL *y, REAL *f^) { f^[0] = y^[1]; f^[1] = -y^[0] * y^[0] * y^[0]; } /*********************************************************************** * alpha-numeric description of^ above system * ***********************************************************************/ static char *dgltxt10(void) { return "y1' = y2\n" "y2' = -y1^3\n"; } /* ----------------------- DE system number 11 ---------------------- */ /* --------------------- (stiff IVP for gear()) --------------------- */ /*********************************************************************** * Compute the value f of the right hand side of a DE system at (x,y). * ***********************************************************************/ static void dgl11(REAL x, REAL *y, REAL *f^) { #def^ine dk ((REAL)-10000.0) f^[0] = y^[1]; f^[1] = dk * y^[0] + (dk - ONE) * y^[1]; #undef^ dk } /*********************************************************************** * alpha-numeric description of above system * ***********************************************************************/ static char *dgltxt11(void) { return "y1' = y1 (dk = -10000)\n" "y2' = dk * y1 + (dk - 1) * y2\n"; } /*********************************************************************** * Compute the value of the analytic solution y(x) for the above DE * * for the initial value problem (y0) = (1,1) * ***********************************************************************/ static void loesung11(REAL x, REAL *y) { y^[0] = EXP(-x); y^[1] = -EXP(-x); } {end of file t_dgls.pas}