{***********************************************************************
*                                                                      *
* Solve a two point boundary problem of first order with the shooting  *
* -------------------------------------------------------------------  *
* method                                                               *
* ------                                                               *
*                                                                      *
* Programming language: ANSI C                                         *
* Author:               Klaus Niederdrenk (FORTRAN)                    *
* Adaptation:           Juergen Dietel, Computer Center, RWTH Aachen   *
* Source:               existing C, Pascal, QuickBASIC and FORTRAN     *
*                       codes                                          *
* Date:                 6.2.1992; 10.30.1995                           *
*                                                                      *
*                                                                      *
*                                  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" [BIBLI 11].           *
***********************************************************************}
UNIT Urwp;

INTERFACE

Uses WinCrt1, Fgauss, Uawp, Ab_mou, Bulirsch;

Type ar_real = DOUBLE;

Function rwp  {Shooting method for boundary value problem of 1st order}
       (
         a,                      {left end point .....................}
         b,                      {right end point ....................}
         h: ar_real;             {starting step size .................}
         y_start: pVec;          {initial approximation or solution ..}
                                 {initial value problem  y(a) ........}
         n: integer;             {number of differntial equations ....}
         awpnumm: integer;       {Number of desired IVP solver .......}
         epsawp,                 {error bound for initial value       }
                                 {problem ............................}
         epsrb: ar_real;         {error bound for boundary value      }
                                 {problem ............................}
         fmax,                   {maximal number of calls of dgl() ...}
         itmax: integer;         {maximal number of Newton iterations }
     var act_iter: integer       {actual number of Newton steps ......}
       ): integer;               {error code .........................}


IMPLEMENTATION

{print a vector of length n and its name}
Procedure zeig(s:String; v:pVec; n:integer);
Var i:integer;
Begin
  writeln(s);
  for i := 0 to n-1 do write(' ', v^[i]);
  writeln
End;

Procedure randbed(ya:pVec; yb:pVec; r:pVec);
Begin
  r^[0] := ya^[0];
  r^[1] := yb^[0]
End;

Function Power(y,x:ar_real):ar_real;
Begin
  Power:=0.0;
  if y<0.0 then exit
  else Power:=Exp(x*Ln(y))
End;

Procedure Copy_Vector(a,b:pVec;n:integer);
Var i:integer;
Begin
  For i:=0 to n-1 do a^[i]:=b^[i]
End;

{ ------------------------------------------------------------------- }
Function rwp  {Shooting method for boundary value problem of 1st order}
       (
         a,                      {left end point .....................}
         b,                      {right end point ....................}
         h: ar_real;             {starting step size .................}
         y_start: pVec;          {initial approximation or solution ..}
                                 {initial value problem  y(a) ........}
         n: integer;             {number of differntial equations ....}
         awpnumm: integer;       {Number of desired IVP solver .......}
         epsawp,                 {error bound for initial value       }
                                 {problem ............................}
         epsrb: ar_real;         {error bound for boundary value      }
                                 {problem ............................}
         fmax,                   {maximal number of calls of dgl() ...}
         itmax: integer;         {maximal number of Newton iterations }
     var act_iter: integer       {actual number of Newton steps ......}
       ): integer;               {error code .........................}

{***********************************************************************
* This function solves the first order boundary value problem          *
*                                                                      *
*     y' = F(x,y),    a <= x <= b,    R(y(a),y(b)) = 0.                *
*                                                                      *
* It uses the shooting method starting with an approximate y_start for *
* the associated initial value y(a), from which it constructs an       *
* approximate solution using as initial value solver the function awp()*
* The nonlinear system arising in the shooting method is solved using  *
* Newton's method.                                                     *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* a        left end point                                              *
* b        right end point  (b > a)                                    *
* h        suitable step size for solving the associated initial value *
*          problem for the shooting method approximately               *
* y_start  [0..n-1] vector, initial approximation for y(a)             *
* n        number of differential equations                            *
* dgl      pointer to the function which evaluates the right hand side *
*          of the differential equation  y' = f(x,y)                   *
* rand     pointer to the function that evaluates the boundary cond.   *
* awpnumm  label for the desired IVP solver in the shooting method:    *
*          = 1: Runge-Kutta embedding formula of 4/5th order; England  *
*               formula from awp().                                    *
*          = 2: Predictor-corrector method of order 4 by Adams-        *
*               Bashforth-Moulton (from prae_korr())                   *
*          = 3: Runge-Kutta embedding formula of 7/8th order (from     *
*               einb_rk())                                             *
*          = 4: Extrapolation method of  Bulirsch-Stoer (from          *
*               bul_stoe())                                            *
*          = 5: implicit Runge-Kutta-Gauss method (from implruku())    *
* epsawp   desired accuracy for the solution of the associated initial *
*          value problem                                               *
* epsrb    desired accuracy for which the approximation y_start for    *
*          the initial value  y(a) of a solution of the boundary value *
*          problem should satisfy the boundary condition R.            *
* fmax     upper bound for number of calls of dgl() in the system of   *
*          differential equations when solving the associated initial  *
*          value problem                                               *
* itmax    upper bound for number of Newton iterations when solving    *
*          the nonlinear systems in the shooting method                *
*                                                                      *
* Output parameters:                                                   *
* ==================                                                   *
* y_start   [0..n-1] approximation for the initial value y(a) of a     *
*           solution  y of the boundary problem                        *
* act_iter  number of Newton iterations performed                      *
*                                                                      *
* Return value :                                                       *
* ==============                                                       *
* error code                                                           *
* = 0: all is ok                                                       *
* = 1: at least one of the accuracy bounds  epsawp, epsrb too small    *
* = 2: b <= a                                                          *
* = 3: h <= 0                                                          *
* = 4: n <= 0                                                          *
* = 5: improper input for awpnumm                                      *
* = 6: maximal number of allowed function evaluations exceeded         *
* = 7: act_iter > itmax. number of allowd Newton steps exceeded without*
*      finding a suitable value for  y_start.                          *
* = 8: The Jacobi matrix for Newton iterations is singular.            *
* = 9: lack of memory space                                            *
* >  9: error in one of the IVP solvers at first node                  *
* > 19: error in one of the IVP solvers at second node                 *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* ar_real, dglsysfnk, rndbedfnk, max, MACH_EPS, copy_vector, FABS,     *
* vminit, vmalloc, vmcomplete, vmfree, VEKTOR, VVEKTOR, MATRIX, awp,   *
* gauss, POW, ZERO, ONE, HALF, prae_korr, einb_rk, bul_stoe, implruku  *
***********************************************************************}
Label fin;
Const
      SINGU=1;             { Name for the return value of gauss(),    }
                           { which indicates that the matrix is       }
                           { singular                                 }
      ENGL45=6;            { awp() shall use the England embedding    }
                           { formula of 4th and 5th order             }

      FORMEL11=11;         { Macros for calling einb_rk(); to improve }
      MIT_HULL=1;          { readability                              }
      VON_VORNE=TRUE;
      NOT_SAVE=FALSE;

      MMAX=12;             { for calling implruku()                   }
      STUETZ='stuetz';
      AUS='';
      PROT='';

      HALF = 0.5;

Var
        Debug:Boolean;     { if True, print debugging infos            }
        yk,                { [0..n-1] vector with value of solution    }
                           { at right end point                        }
        yaj,               { [0..n-1] vector of modified y_start used  }
                           { for forming the Jacobi matrix             }
        r,                 { [0..n-1] vector with left boundary        }
                           { condition, which serves as the right hand }
                           { side for the linear system in the Newton  }
                           { method.                                   }
        rj,                { [0..n-1] vector with left end boundary    }
                           { condition of the modified initial value   }
                           { problem.                                  }
        d: pVec;           { [0..n-1] vector, the solution of the      }
                           { linear system.                            }
        amat: pMat;        { [0..n-1,0..n-1] array of the Jacobi       }
                           { matrix for the  Newton step.              }
        yk2,               { [0..n-1] vector with the values of the    }
                           { solution of the DE at the left endpoint   }
                           { for implruku().                           }
        g: pVec;           { [0..n-1] vector of weights for implruku() }
        xk,                { left end point for  awp()                 }
        hk,                { desired step size for awp()               }
        epsabs,            { absolute accuracy desired in  awp()       }
        epsrel,            { relative accuracy in awp()                }
        epsrel2,           { relative error bound for implruku()       }
        epsrel3,           { same as `epsrel2', but possibly altered   }
                           { during last call of  implruku()           }
        rmax,              { Maximum norm of the left boundary         }
                           { condition                                 }
        delta,             { Aux variable for  Jacobi matrix           }
        mach1,             { accuracy bounds depending on computer     }
        mach2: ar_real;    {                                           }
        aufrufe:integer;   { number of function evaluations in awp()   }
        pivot: pIVec;      { Permutation vector for Gaussian elimin.   }
        i,                 { loop counter                              }
        jacobi,            { counter for columns of Jacobi matrix      }
        mark,              { error code from gauss()                   }
        sign_det,          { sign of the permutation in gauss()        }
        fehler: integer;   { error code in  awp()                      }

Begin

  { --------- allocate storage for aux vectors ------------------------- }
  New(yk); New(yaj); New(r); New(rj); New(d); New(amat); New(yk2); New(g);
  New(pivot);

  mach1  := Power(MACH_EPS, 0.75);
  mach2  := 100.0 * MACH_EPS;

  Debug:=False;

  if epsrb < mach2 then rwp:=1;               {check input}
  if b <= a then rwp:=2;
  if h <= mach2 * ABS(a) then rwp:=3;
  if n <= 0 then rwp:=4;
  if (awpnumm < 1) or (awpnumm > 5) then rwp:=5;

  act_iter := 0;
  epsabs    := HALF * epsawp;
  epsrel    := epsabs;

  for i := 0 to n-1 do g^[i]:=ONE;   {weight vector for implruku(), set}
                                     {up as all ones.                  }

  if epsrel < 1E-10 then epsrel2:=1e-10 else epsrel2 := epsrel;

  {*********************************************************************
  * If  y_start proves to be a sufficiently good approximation for y(a)*
  * in the following loop, we report the success and check the         *
  * remaining input data in awp().                                     *
  *********************************************************************}

  While True do
  begin

    if Debug then zeig('Y_start:', y_start, n);

    copy_vector(yk, y_start, n);
    xk := a;
    hk := h;

    Case awpnumm of
      1: awp(xk, b, n, yk, epsabs, epsrel, hk, ENGL45, fmax, aufrufe, fehler);
      2: fehler := prae_korr(xk, yk, n, b, hk, epsabs,epsrel, fmax, aufrufe, ONE,VON_VORNE);
      3: fehler := bul_stoe(xk, b, n, yk, epsabs, epsrel,hk, ONE, VON_VORNE, fmax, aufrufe)
    end;

    if fehler <> 0 then
    begin
      rwp := 10+fehler;      {error in IVP solver ? }
      goto fin
    end;

    if Debug then zeig('Yk:', yk, n);

    randbed(y_start, yk, r);

    rmax:=0.0;
    for i := 0 to n-1 do rmax := Max(rmax, ABS(r^[i]));
    if rmax < epsrb then              {boundary condition satisfied ? }
    begin                                     {report success         }
      rwp := 0;
      goto fin
    end;

    Inc(act_iter);
    if act_iter > itmax then      {approximation not good enough      }
    begin                         {after  itmax Newton steps ?        }
      rwp := 6;                   {report failure                     }
      goto fin
    end;

    {find a better approximation  y_start for y(a) by using the
     Newton method with a Jacobi matrix that is approximated by one-
     sides difference quotients    }

    for jacobi := 0 to n-1 do
    begin
      copy_vector(yk,  y_start, n);
      copy_vector(yaj, y_start, n);

      if ABS(yk^[jacobi]) < mach2 then
      begin
        yk^[jacobi] := yk^[jacobi] + mach1;
        delta       := ONE / mach1
      end
      else
      begin
        yk^[jacobi] := yk^[jacobi] * (ONE + mach1);
        delta       := ONE / (mach1 * yk^[jacobi])
      end;

      yaj^[jacobi] := yk^[jacobi];
      xk := a;
      hk := h;

      if Debug then zeig('Yaj:', yaj, n);

      Case awpnumm of
        1: awp(xk, b, n, yk, epsabs, epsrel, hk, ENGL45, fmax, aufrufe, fehler);
        2: fehler := prae_korr(xk, yk, n, b, hk, epsabs, epsrel, fmax, aufrufe, ONE, VON_VORNE);
        3: fehler := bul_stoe(xk, b, n, yk, epsabs, epsrel, hk, ONE, VON_VORNE, fmax, aufrufe)
      End;

      if fehler <> 0 then                      {error in IVP solver ? }
      begin
        rwp := 20 + fehler;
        goto fin
      end;

      if Debug then zeig('Yk:', yk, n);

      randbed(yaj, yk, rj);

      for i := 0 to n-1 do
        amat^[i][jacobi] := (rj^[i] - r^[i]) * delta

    end;

    gauss(0, n, amat, amat, pivot, r, d, sign_det, mark);

    if mark = SINGU then                   {Jacobi matrix singular? }
    begin                                  {return error            }
      rwp := 8;
      goto fin
    end;                             

    if Debug then zeig('D:', d, n);

    for i := 0 to n-1 do                           {correct y_start }
      y_start^[i] := y_start^[i] - d^[i];

  end; {while True}
  {free memory}
fin:Dispose(yk); Dispose(yaj); Dispose(r); Dispose(rj); Dispose(d); Dispose(amat);
    Dispose(yk2); Dispose(g); Dispose(pivot);
End;

END.

{ ---------------------------- END rwp.pas ----------------------------}