###############################################################

##

## Newton's method for Maple

##

##    Procedure for solving systems of nonlinear equations

##    using Newton's method.

##

##    Written by Ross Taylor, Clarkson University, Potsdam, New York.

##    November-December 1993

##    Revised June 1994

##    Revised June 1995 to include timelimit option

##    Revised February 1997 to deal with functions as procedures

##    Revised March 1998 by Richard Baur to deal with hardware floating point computation

##

`help/text/Newton` := TEXT (

`FUNCTION: Newton - iterative method for obtaining numerical `,

`                   solutions to systems of nonlinear equations.`,

``,

`CALLING SEQUENCE:`,

`   Newton (Eqns, Varlist, <options>);`,

`  `,

`PARAMETERS:`,

`   Eqns    -   Vector or list of equations (or a single equation)`,

`   Varlist -   Vector or list of equations (or just one equation)`,

`               of the form name = numeric, where name is the name`,

`               of an unknown variables and numeric is the initial`,

`               estimate of that variable.`,

``,

`SYNOPSIS:`,

`- Newton is a fairly straightforward implementation of Newton's`,

`  method for finding numerical solutions to systems of nonlinear`,

`  algebraic equations.  Newton is somewhat more flexible than Maple's`,

`  built-in solver, fsolve, by permitting users to control the`,

`  iteration history by (i) providing an initial estimate of all`,

`  variables (required), and (ii) providing a procedure that `,

`  controls the change in variables between iterations (optional).`,

``,

`- Optional arguments must be given in the form`,

`    option = something`,

`  Valid options are: iterations, output, tolerance, timelimit, and`,

`  steps and are discussed below.`,

``,

`- If a converged solution is not found within a specified number`,

`  of iterations, Newton will fail with an error message to that`,

`  effect.  The default number of iterations is 25 but that can be`,

`  changed by the optional argument iterations = x, where x is an`,

`  integer greater than 1.`,

``,

`- If a converged solution is not found within a specified time,`,

`  Newton will fail with an error message to that effect.  The`,

`  default time limit is 3600 seconds but that can be changed`,

`  by the optional argument timelimit = x, where x is a number`,

`  of seconds. If you would like Newton to run indefinitely, then`,

`  use timelimit = infinity.`,

``,

`- Newton can provide a complete iteration history.  This can be `,

`  useful in diagnosing problems with examples that fail to converge.`,

`  By default Newton does not provide any intermediate output which`,

`  is controlled by the output = {} option where {} is a set (or list)`,

`  of keywords.  Allowed keywords are:`,

`     norm      - prints iteration count and function norm on each iteration`,

`     variables - prints variable values on each iteration`,

`     functions - prints function values`,

`     jacobian  - prints jacobian matrix`,

`     sparsity  - prints sparsity pattern of jacobian`,

`     symbolics - prints equations and jacobian element expressions`,

``,

`- Newton will perform iterations until the 2-norm (use help for norm if you`,

`  you don't know what that is) is less than some small number.  The value`,

`  assigned by Newton is 10^(2-Digits).  The tolerance can be changed using`,

`  the option tolerance = x, where x is a small number. Newton prints a `,

`  warning message if the tolerance is smaller than the default value. To `,

`  converge to a higher accuracy it is necessary to make Digits larger.`,

``,

`- Newton's method will not always converge.  Some problems are very hard`,

`  to solve even from very good initial estimates.  In some cases Newton's`,

`  method exhibits periodic cycles or even chaotic iteration histories.`,

`  It will not always be possible to solve such problems with Newton's`,

`  method but it is sometimes possible to alleviate convergence problems`,

`  by controlling the Newton-step on each iteration.  Examples of step`,

`  controls that may be useful include preventing certain variables from`,

`  taking on meaningless values (for example, many parameters in problems`,

`  with a physical origin cannot be negative).  Alternatively, it might be`,

`  desirable to limit the Newton-step to some maximum value.  Newton`,

`  can use a user-supplied procedure for controlling the Newton steps.`,

`  To indicate that a procedure is to be used it is necessary to use the`,

`  option steps = procname, where procname is the name of the procedure.`,

`  The procedure may contain up to three arguments in order:`,

`      vector of original step changes computed in Newton`,

`      vector of current variable values`,

`      vector of variable names.`,

`  See the section entitled EXAMPLES for illustrations.`,

``,

``,

`EXAMPLES:`,

``,

`> with(linalg):`,

`> Newton(3*x-exp(-x)=0,x=1,output={variables,functions});`,

`            x = 1`,

`            f[1] = 3 - exp(-1)`,

`            x = .2184635450`,

`            f[1] = -.1483621447`,

`            x = .2574676952`,

`            f[1] = -.0006035121`,

`            x = .2574676504`,

``,

`> f1 := x1^2 + x2^2 = 17:`,

`> f2 := (8*x1)^(1/3) + sqrt(x2) = 4:`,

`> Newton ([f1,f2], [x1=0.8, x2=5.8]);`,

`            [x1=0.99999999998, x2=4.0000000000]`,

``,

`> Newton ([f1,f2], [x1=0.8, x2=5.8],timelimit=0.1);`,

`  Error, (in Newton) Time limit exceeded`,

``,

`Here we give examples of procedures for controlling the Newton step.  `,

`The procs shown below are part of the Newton package and were loaded`,

`with Newton.`,

``,

`The first is a proc for simply damping (or accelerating) the Newton step.`,

`Nothing fancy like line searching here.  Just multiplication of the`,

`Newton step vector with a user determined (global) number.`,

``,

`> dampstep := proc(delx)`,

`  # Procedure for limiting Newton steps (called from Newton)`,

`  # delx - vector of original steps`,

`  # dampfactor is a global vriable set outside the proc so that user `,

`  # need not pass a value by call to Newton.`,

`  local i, newstep, nvars;`,

`  if not type(dampfactor,numeric) then `,

`    ERROR(``Please set global variable named dampfactor before using dampstep``);`,

`  fi;`,

`  nvars := linalg[vectdim](delx);`,

`  newstep := linalg[vector](nvars);`,

`  for i from 1 to nvars do `,

`    newstep[i] := dampfactor * delx[i];`,

`  od;`,

`  evalm(newstep);`,

`  end;`,

``,

`The following proc limits the correction to some maximum step size that the`,

`user determines.`,

``,

`> maxchange := proc(delx, dxmax)`,

`  # Procedure to limit the maximum step size`,

`  # delx - original step`,

`  # dxmax - maximum step size (must be positive number)`,

`  if delx > dxmax then`,

`    dxmax `,

`  elif delx < -dxmax then`,

`    -dxmax `,

`  else `,

`    delx `,

`  fi; `,

`  end:`,

``,

`Now a routine that uses maxchange.`,

``,

`> teststep := proc(delx)`,

`  # Procedure for limiting Newton steps`,

`  # delx - vector of original steps`,

`  # maxstep is a global variable set outside the proc `,

`  # so that user need not pass a value by call to Newton`,

`  local i, nvars, newstep;`,

`  if not type(maxstep,numeric) then `,

`     ERROR(``Please set global variable named maxstep before using teststep``); `,

`  fi;`,

`  nvars := linalg[vectdim](delx);`,

`  newstep := linalg[vector](nvars);`,

`  for i from 1 to nvars do `,

`    newstep[i] := maxchange(delx[i], maxstep); `,

`  od;`,

`  evalm(newstep);`,

`  end:`,

``,

`> maxstep := 0.1:`,

`> Newton ([f1,f2], [x1=0.8, x2=5.8],steps=teststep);`,

`            [x1=4.071504896, x2=.6502675557]`,

``,

`Note that Newton converged to a different result this time.`,

``,

`The next three procedures prevent variables from going outside some user`,

`defined region.  In the first of these procs, the variable is constrained `,

`between a lower and an upper limit.  In the subsequent examples only one`,

`of these limits is imposed.  These procs would be called from a proc `,

`similar to teststep above.`,

``,

`> trustregion := proc (x, delx, xmin, xmax) `,

`  local temp; `,

`  temp := x-delx; `,

`  if temp < xmin then `,

`    RETURN (1/2*x-1/2*xmin);`,

`  elif xmax < temp then `,

`    RETURN (1/2*x-1/2*xmax); `,

`  else `,

`    RETURN (delx);`,

`  fi; `,

`  end;`,

` `,

`> lowerbound := proc (x, delx, xmin) `,

`  local temp; `,

`  temp := x-delx; `,

`  if temp < xmin then `,

`    RETURN (1/2*x-1/2*xmin);`,

`  else `,

`    RETURN (delx);`,

`  fi; `,

`  end;`,

`  `,

`> upperbound := proc (x, delx, xmax) `,

`  local temp; `,

`  temp := x-delx; `,

`  if xmax < temp then `,

`    RETURN (1/2*x+1/2*xmin);`,

`  else `,

`    RETURN (delx);`,

`  fi; `,

`  end;`,

``

):













###############################################################

##

##  Main procedure

##



Newton := proc (Eqns1:: {equation,list,vector,procedure},

                 Vars:: {equation,list(name=complex)})



local eps,          # tolerance

      infolevel,    # used to control printed output

      funcJ,        # array containing jacobian element functions

      maxit,        # maximum number of iterations

      ProcPath,     # functionname prefix

      nextstep,     # name of proc used to control steps

      a,s,t,        # Maximum allowed run time (seconds)

      Otherargs,    # List of arguments to pass to function procedure

      ShowTime,     # Switch for printing the time at the end

      comTime,      # Computation time

      res,          # result

      k,            # help variable

      maxtime:      # Maximum allowed run time (seconds)





####

# look over optional arguments and set parameters

eps := 10^(-Digits+2);

infolevel := array(1..7, [0,0,0,0,0,0,0]);

maxit := 25;  # Default number of iterations

maxtime := 3600; # Default time limit in seconds

nextstep := SimpleStep; # simple procedure not modifying the change of the

                        #  the state vector

Otherargs := 0;

ShowTime := false:

ProcPath := `NEWTON/`:

for a in [args[3..nargs]] do

        if not type(a,equation) then

          ERROR(`Incorrect optional argument`,a)

        fi;

        s := lhs(a); t := rhs(a);

        if s = 'tolerance'  then

          eps := t;

          if eps < 10^(-Digits+2) then

            lprint(`WARNING: tolerance may be too small`)

          fi;

        elif s = 'output' then

            if type(t,list) or type(t,set) then

                if member(`norm`,t) then infolevel[1] := 1 fi;

                if member(`variables`,t) then infolevel[2] := 1 fi;

                if member(`functions`,t) then infolevel[3] := 1 fi;

                if member(`jacobian`,t) then infolevel[4] := 1 fi;

                if member(`sparsity`,t) then infolevel[5] := 1 fi;

                if member(`symbolics`,t) then infolevel[6] := 1 fi;

                if member(`time`,t) then ShowTime := true fi;

            fi;

            for k in select(type, t ,`equation`) do

             if lhs(k) = `variables` then

                 infolevel[2] := -1:

                 infolevel[7] := rhs(k):

              fi:

            od:

        elif s = 'iterations' then

            maxit := t;

        elif s = 'timelimit' then

            maxtime := t;

        elif s = 'steps' then

            nextstep := t;

        elif s = jacobian then

            if type(t,procedure) then

              funcJ := t;

            elif t = numeric then

              funcJ := 0;

            elif t = auto then

              lprint(`WARNING:: not implemented yet`):

            fi;

        elif s = otherargs then

            Otherargs := op(t);

        elif s = mode then

            if t = hf then

              ProcPath := `NEWTON/hf/`:

            fi:

        fi;

od;





if ShowTime then comTime := time(): fi:



if ProcPath  = `NEWTON/hf/` then

 res := `NEWTON/hf/Newton`(Eqns1, Vars, funcJ, eps, infolevel,

                              maxit, nextstep, maxtime, ProcPath,

                              Otherargs);

else



 res := `NEWTON/Newton`(Eqns1, Vars, funcJ, eps, infolevel,

                              maxit, nextstep, maxtime,

                              Otherargs);

fi:



if ShowTime then

 print(cat(` Time needed for the computation:`,

            convert(time()- comTime, string), ` sec`)):

fi:



RETURN(res):



end;

###############################################################

###############################################################



###############################################################

##

## Newton routine based on the previous version

##



`NEWTON/Newton` := proc (Eqns1, Vars, funcJ, eps, infolevel,

                              maxit, NextStepProc, maxtime,

                              Otherargs)



local neqns,        # number of equations

      nvars,        # number of variables

      Fvec,         # vector of functions

      Ftemp,        # temporary vector of functions

      F1,           # vector of functions used for numerical Jacobian

      Jmat,         # array containing jacobian element functions

      i, j, k,      # various counters

      dfij,         # used to flag zero in jacobian matrix

      delx,         # change in x vector

      del1,         # adjusted change in x vector

      nrm,          # norm of function vector

      xnrm,         # norm of variable change vector

      nextstep,     # name of proc used to control steps

      AddArgs  ,    # List of arguments to pass to function procedure

      ff,           # function vector

      JJ,           # jacobian matrix

      itercount,    # records current iteration number

      a,s,t,        # used in processing optional arguments

      x,            # vector of variable names

      xx,           # vector of variable values

      xstart,       # vector of initial values for entries in x

      xstep,        # step size for numerical Jacobian

      Eqns,         # copy of Eqns1 in argument list

      Varlist,      # copy of Vars in argument list

      Varseq,       # Sequence of variables

      leftovernames,# Unitialised variables

      timenow,      # Current time (s)

      starttime,    # Time the calculations started (s)

      et;           # Elapsed time







## Some of the arguments have to be redefined:

if Otherargs = 0 then

 AddArgs := NULL:

else

 AddArgs := Otherargs:

fi:

if NextStepProc = SimpleStep then

 nextstep := NULL:

else

  nextstep := NextStepProc:

fi:



# Determine number of equations

if type (Eqns1, equation) then

  Eqns := [Eqns1];

elif type (Eqns1,procedure) then

else

  Eqns := Eqns1;

fi;



if type (Vars, equation) then

  Varlist := [Vars];

else

  Varlist := Vars;

fi;



#neqns := vectdim(Eqns);

nvars := linalg[vectdim](Varlist);



#if neqns <> nvars then

#   ERROR(`Numbers of equations and variables are not the same`);

#fi;



x := linalg[vector](nvars);

xstart := linalg[vector](nvars);



for i from 1 to nvars do

     x[i] := lhs(Varlist[i]);

     xstart[i] := rhs(Varlist[i]);

od;



# Define function vector and jacobian matrix of appropriate size

Fvec:=linalg[vector](nvars):

Jmat := array(sparse,1..nvars,1..nvars):



# Create variable sequence

Varseq := seq(x[k],k=1..nvars);



if not type(Eqns1,procedure) then

  leftovernames := indets({op(Eqns)},name) minus {Varseq};

  if leftovernames <> {} then

    ERROR(`The following variable(s) is/are not initialized `,leftovernames);

  fi;

fi;



# Create function vector

if not type (Eqns1,procedure) then

  for i from 1 to nvars do

    if type(Eqns[i],`=`) then

      Ftemp:=(lhs(Eqns[i])-rhs(Eqns[i]));

    else

      Ftemp:=(Eqns[i]);

    fi;

    Fvec[i]:=unapply(Ftemp,Varseq);

    for j from 1 to nvars do

      if has(Ftemp,x[j]) then

        dfij := diff(Ftemp,x[j]);

        Jmat[i,j] := unapply(dfij,Varseq);

      fi;

    od;

  od;

else

#   print(`Jmat already assigned`);

  Fvec:=Eqns1;

fi;





if infolevel[6] > 0 then print(Fvec(Varseq,AddArgs)); fi;



if infolevel[6] > 0 then print(Jmat(Varseq)); fi;

if infolevel[5] > 0 then

  print(plots[sparsematrixplot](Jmat(Varseq),title=`Sparsity pattern of jacobian (upside down)`));

fi;



# Define a few vectors and matrices for working space

xx := linalg[vector](nvars):

xx := xstart;

delx := linalg[vector](nvars);

ff := linalg[vector](nvars);

F1 := linalg[vector](nvars);

del1 := linalg[vector](nvars);

JJ := array(sparse,1..nvars,1..nvars);



# Evaluate the function vector and jacobian matrix and the norm of the function vector

ff:=Fvec(seq(xx[k],k=1..nvars),AddArgs);

if type(Jmat,procedure) then

   JJ:=Jmat(seq(xx[k],k=1..nvars),AddArgs);

elif type(Eqns1,procedure) and not type(Jmat,procedure) then

  for j from 1 to nvars do

      xstep := max(1e-4*abs(xx[j]),1e-10);

      xx[j] := xx[j] + xstep;

      F1 := Fvec(seq(xx[k],k=1..nvars),AddArgs);

      for i from 1 to nvars do

          JJ[i,j] := (F1[i] - ff[i]) / xstep;

      od;

      xx[j] := xx[j] - xstep;

  od;

else

   JJ:=Jmat(seq(xx[k],k=1..nvars));

fi;



nrm := evalf(linalg[norm](ff,2));

xnrm := nrm;





# Perform iterations

itercount := 0;

starttime := time();

while (nrm > eps and xnrm > eps) do



   # Print out various things according to inflolevel

   if infolevel[1] > 0 then print(`Iteration `. itercount, `   Norm = `, nrm, `   Xnrm = `,xnrm); fi;

   if infolevel[2] > 0 then print(seq(x[k]=xx[k],k=1..nvars)); fi;



   if infolevel[2] < 0 then print(

                        select(has,[seq(x[k]=xx[k],k=1..nvars)],

                               infolevel[7]));

   fi;

   if infolevel[3] > 0 then print(seq(`f`[k]=ff[k],k=1..nvars)); fi;

   if infolevel[4] > 0 then print(JJ); fi;



   # Solve linear system to get vector of corrections

   delx:=linalg[linsolve](JJ,ff);



   # Modify steps if required by optional arguments

   if nextstep <> NULL then

      del1:=nextstep(delx,xx,x);

   else

      del1 := delx;

   fi;



   # Compute new set of variable values and the next function vector etc.

   xx:=evalf(evalm(xx-del1));

   ff:=Fvec(seq(xx[k],k=1..nvars),AddArgs);

   if type(Jmat,procedure) then

      JJ:=Jmat(seq(xx[k],k=1..nvars),AddArgs);

   elif type(Eqns1,procedure) and not type(Jmat,procedure) then

     for j from 1 to nvars do

         xstep := max(1e-6*abs(xx[j]),1e-10);

         xx[j] := xx[j] + xstep;

         F1 := Fvec(seq(xx[k],k=1..nvars),AddArgs);

         for i from 1 to nvars do

           if abs(F1[i] - ff[i]) < 1e-20 then

            JJ[i,j] := 0:

           else

            JJ[i,j] := (F1[i] - ff[i]) / xstep;

           fi:

         od;

         xx[j] := xx[j] - xstep;

     od;

   else

      JJ:=Jmat(seq(xx[k],k=1..nvars));

   fi;

#   JJ:=Jmat(seq(xx[k],k=1..nvars));

   nrm := evalf(linalg[norm](ff,2));

   xnrm := evalf(linalg[norm](delx,2));



   # Increment iteration count and check that we have not done too many iterations

   itercount := itercount+1;

   if itercount > (maxit + 1) then ERROR(`Maximum number of iterations exceeded`);  fi;



   # Check time

   timenow := time();

   et := timenow - starttime;

   if et > maxtime then ERROR(`Time limit exceeded`); fi;

od:



# If we have made it this far we probably have a converged solution

if nvars=1 then

  RETURN(x[1]=xx[1]);

else

  RETURN([seq(x[k]=xx[k],k=1..nvars)]);

fi;

end:



#

# Here we give some procedures for controlling the Newton step.

# The first limits the correction to some maximum step size that

# the user determines.  The use of these routines will become

# clear later.



maxchange := proc(delx, dxmax)

# Procedure to limit the maximum step size

# delx - original step

# dxmax - maximum step size (must be positive number)

if delx > dxmax then

  dxmax elif delx < -dxmax then-dxmax

else

  delx

fi;

end:



# Now a routine that uses maxchange.



teststep := proc(delx)

# Procedure for limiting Newton steps (called from Newton)

# delx - vector of original steps

#           proc returns a vector of steps modified accoring to the rules in this proc

# maxstep is a global vriable set outside the proc so that user need not pass a value by call to Newton

local i, nvars, newstep;

if not type(maxstep,numeric)  then ERROR(`Please set global variable named maxstep before using teststep`); fi;

nvars := linalg[vectdim](delx);

newstep := linalg[vector](nvars);

for i from 1 to nvars do newstep[i] := maxchange(delx[i], maxstep); od;

evalm(newstep);

end:



# Here is a proc for simply damping (or accelerating) the Newton step.

# Nothing fancy like line searching here.  Just multiplication of the

# Newton step vector with a user determined (global) number.



dampstep := proc(delx)

# Procedure for limiting Newton steps (called from Newton)

# delx - vector of original steps

# dampfactor is a global vriable set outside the proc so that user

# need not pass a value by call to Newton.

local i, newstep, nvars;

if not type(dampfactor,numeric) then ERROR(`Please set global variable named dampfactor before using dampstep`); fi;

nvars := linalg[vectdim](delx);

newstep := linalg[vector](nvars);

for i from 1 to nvars do newstep[i] := dampfactor * delx[i]; od;

evalm(newstep);

end:



trustregion := proc (x, delx, xmin, xmax)

local temp;

temp := x-delx;

if temp < xmin then

  RETURN (1/2*x-1/2*xmin);

elif xmax < temp then

  RETURN (1/2*x-1/2*xmax);

else

  RETURN (delx);

fi;

end:



lowerbound := proc (x, delx, xmin)

local temp;

temp := x-delx;

if temp < xmin then

  RETURN (1/2*x-1/2*xmin);

else

  RETURN (delx);

fi;

end:



upperbound := proc (x, delx, xmax)

local temp;

temp := x-delx;

if xmax < temp then

  RETURN (1/2*x+1/2*xmin);

else

  RETURN (delx);

fi;

end:



###############################################################

##

## Procedure for hardware floating point calculations

##

`NEWTON/hf/Newton` := proc(Eqns1, Vars, funcJ, eps, infolevel,

                              maxit, nextstep, maxtime, ProcPath,

                              Otherargs)

local neqns,        # number of equations

      nvars,        # number of variables

      Fvec,         # vector of functions

      Jmat,         # array containing jacobian element functions

      i, j, k,      # various counters

      dfij,         # used to flag zero in jacobian matrix

      a,s,t,        # used in processing optional arguments

      x,            # vector of variable names

      xstart,       # vector of initial values for entries in x

      Eqns,         # copy of Eqns1 in argument list

      Varlist,      # copy of Vars in argument list

      Varseq,       # Sequence of variables

      InternalSeq,  # Internal sequence of variables

      leftovernames,# Unitialised variables

      DUMMYeq,      # Dummy name

      StepProc,     # Next step procedure

      VarPrint,     # Help variable to store the indices of the

                    # variables for printing information

      VarTemp,      # Temporary help variable

      funcF:        # Procedure assigning the function vector



####

# Some definitions required for the iteration



# Modifications for a single equation

if type (Eqns1, equation) then

  Eqns := [Eqns1];

else

  Eqns := Eqns1;

fi;

if type (Vars, equation) then

  Varlist := [Vars];

else

  Varlist := Vars;

fi;



nvars := nops(Varlist);



if infolevel[2] < 0 then

 VarPrint := array(1..nops(infolevel[7])):

 infolevel[7] := [op(infolevel[7])]:

 VarTemp :=  subs(seq(lhs(Varlist[ii]) = ii, ii = 1..nvars),

                  infolevel[7] ):

 lprint(`Legend:`):

 for k from 1 to  nops(infolevel[7]) do

  if type(VarTemp[k],integer) then

   lprint(`  `,

          convert((infolevel[7][k]), string),

          `  is indexed by  `,

          convert(eval(VarTemp[k]), string),

          ` in the state vector`);

   VarPrint[k] := VarTemp[k]

  else

   lprint(`  `,

          convert((infolevel[7][k]), string),

          ` is not a state`);

   VarPrint[k] := -1:

  fi:

 od:

 infolevel[7] := nops(infolevel[7]):

else

  VarPrint := 0;

fi:



#if neqns <> nvars then

#   ERROR(`Numbers of equations and variables are not the same`);

#fi;



x := array(1..nvars);

xstart := array(1..nvars);



for i from 1 to nvars do

     x[i] := lhs(Varlist[i]);

     xstart[i] := rhs(Varlist[i]);

od;



# Define function vector and jacobian matrix of appropriate size

Fvec := array(1..nvars):

Jmat := array(sparse,1..nvars,1..nvars):



# Create variable sequence

Varseq := seq(x[k],k=1..nvars);



####

# Equation system is given as a list

if not type (Eqns,procedure) then

  # Check of the equality of unknowns and equations

  leftovernames := indets({op(Eqns)},name) minus {Varseq};

  if leftovernames <> {} then

    ERROR(`The following variable(s) is/are not initialized `,leftovernames);

  fi;

  # Create internal representation sequence

  InternalSeq :=  seq(x[k] = _AA[k],k=1..nvars);

  # Create function vector and Jacobian matrix

  for i from 1 to nvars do

    if type(Eqns[i],`=`) then

      Fvec[i]:=subs(InternalSeq, lhs(Eqns[i])-rhs(Eqns[i]));

    else

      Fvec[i]:=subs(InternalSeq, Eqns[i]);

    fi;



    # if Jmat is set to 0 then then jacobian will be computed numerical

    #  this will be done in the NewtonShell routine

    # if Jmat is a preocedure then it is predefined by the user

    if (not(funcJ = 0)  and not(type(funcJ, procedure))) then

     for j from 1 to nvars do

       if has(Fvec[i], _AA[j]) then

         Jmat[i,j] := diff(Fvec[i],_AA[j]);

       else

        Jmat[i,j] := 0;

       fi;

     od;

    fi;

  od;



  if infolevel[6] > 0 then print(Fvec(Varseq,Otherargs)); fi;



  # create function for `evalhf` computations:

  funcF :=  readlib(procmake)(

   `&proc`([_AA,DUMMYeq], [], [],

   `&statseq`(seq(

   `&:=`(`&args`[2][i], Fvec[i]),

         i=1..nvars))

   ));



  # if Jmat is set to 0 then then jacobian will be computed numerical

  #  this will be done in the NewtonShell routine

  # if Jmat is a preocedure then it is predefined by the user

 if (not(funcJ = 0)  and not(type(funcJ, procedure))) then

   funcJ := readlib(procmake)(

    `&proc`([_AA,DUMMYeq], [], [],

     `&statseq`(seq(seq(

     `&:=`(`&args`[2][i,j], Jmat[i,j]),

    i=1..nvars),j=1..nvars))

   ));

  fi:



####

# equation system is given in a procedure

else

  #   print(`Jmat already assigned`);

  # limited syntax-check

  `NEWTON/hf/ProcTest/ff`(Eqns, xstart, `TEST_VECTOR`,nvars):



  # definition of the procedure for the function array

  funcF:=Eqns;

  # set the jacobian computation to a finite difference method.

 if not(type(funcJ, procedure)) and (not(funcJ = 0)) then

  funcJ := 0;

 fi:

fi;



####

# limited syntax-check #

if (type(funcJ, procedure)) then

    `NEWTON/hf/ProcTest/ff`(funcJ, xstart, `TEST_MATRIX`, nvars):

fi:



####

# some manipulation for setting the nextstep procedure

#  check for builtin procedures

if member(nextstep,[`SimpleStep`, `maxchange`, `dampstep`]) then

 # the name of the procedure has to be change since they don't

 #  coincide with the ones use in the analytical version of

 #  Newton.



 StepProc := cat(ProcPath, convert(nextstep,string)):

  # some simple checks for global parameters



 if nextstep = `maxchange` then

   if not type(dxmax,numeric) then

      lprint(`WARNING:: >>maxchange<< requires a global limit stored in dxmax `):

      ERROR(`Set global variable named:     dxmax  `);

   fi:

 elif nextstep = `teststep` then

   ERROR(`Procedure teststep not available, please use maxchange`);



 elif nextstep = dampstep then

    if not type(dampfactor,numeric) then

      lprint(`WARNING:: >>dampstep<< requires a global limit stored in dampfactor `):

      ERROR(`Set global variable named:     dampfactor  `);

   fi:

  fi:

else

 if not(type(nextstep, procedure))  then

   ERROR(`NEXTSTEP function`, convert(nextstep, string), ` is not known or defined.`);

 else

   `NEWTON/hf/ProcTest/ns`(nextstep, xstart, nvars):

   StepProc := nextstep:

 fi:

fi;



# limited syntax-check for the step functions



####

# Iteration

evalhf(

 `NEWTON/hf/NewtonShell`

   (funcF, funcJ, var(xstart), nvars, eps, maxtime, maxit,

    infolevel, StepProc, VarPrint)):



####

# Return of the result



if type(Eqns1, equation) then

 RETURN(x[1] = eval(xstart[1])):

else

 RETURN([seq( x[k] = eval(xstart[k]), k = 1..nvars)]):

fi:



######

end:



###############################################################

##  NEWTON step routines

##

## Here we give some procedures for controlling the Newton step.

## The first limits the correction to some maximum step size that

## the user determines.  The use of these routines will become

## clear later.

##

## NOTE the arguments are order in the following way

##

##      proc(<ARRAY OF COMPUTED DEVIATION VALUES>,

##           <ARRAY OF THE ACTUAL STATE

##             WHICH HAS TO BE RETURNED WITH THE NEW VALUES),

##           <NUMBER OF STATE>)

## NOTE the second argument is the array of the actual states, which are

##      modified in the routines!

##



####

# Procedure to limit the maximum step size

# delx - original step

# dxmax - maximum step size (must be positive number)



`NEWTON/hf/SimpleStep` := proc(delx, xx, nvars)

 local k:

 for k to nvars do

  xx[k] := xx[k] - delx[k]:

 od

end:



####

# Procedure to limit the maximum step size

# delx - original step

# dxmax - maximum step size (must be positive number)

`NEWTON/hf/maxchange` := proc(delx, xx, nvars)

local k;

for k to nvars do

 if delx[k] > dxmax then

   xx[k] := xx[k] - dxmax

 elif delx[k] < -dxmax then

   xx[k] := xx[k] + dxmax

 else

   xx[k] := xx[k] - delx[k]

 fi;

od:

end:





# Here is a proc for simply damping (or accelerating) the Newton step.

# Nothing fancy like line searching here.  Just multiplication of the

# Newton step vector with a user determined (global) number.



####

# Procedure for limiting Newton steps (called from Newton)

# delx - vector of original steps

# dampfactor is a global vriable set outside the proc so that user

# need not pass a value by call to Newton.

`NEWTON/hf/dampstep` := proc(delx, xx, nvars)

local k:

for k to nvars do

 xx[k] := xx[k] -  dampfactor * delx[k];

od:

end:



####

####

`NEWTON/hf/trustregion` := proc (x, delx, xmin, xmax)

local temp;

temp := x-delx;

if temp < xmin then

  RETURN (1/2*x-1/2*xmin);

elif xmax < temp then

  RETURN (1/2*x-1/2*xmax);

else

  RETURN (delx);

fi;

end:



`NEWTON/hf/lowerbound` := proc (x, delx, xmin)

local temp;

temp := x-delx;

if temp < xmin then

  RETURN (1/2*x-1/2*xmin);

else

  RETURN (delx);

fi;

end:



`NEWTON/hf/upperbound` := proc (x, delx, xmax)

local temp;

temp := x-delx;

if xmax < temp then

  RETURN (1/2*x+1/2*xmin);

else

  RETURN (delx);

fi;

end:







###############################################################

##

## Procedures for the hf-mode

##

###############################################################



###############################################################

##

##  procedure testing a procedure for evalhf statements

##



`NEWTON/hf/ProcTest/ff` := proc(EqProc, xstart, test, nvars)

  local error_statement,  # error returned by evalhf evaluation

        test_array:       # array for testing the procedure

  if test = `TEST_VECTOR` then

   test_array := array(1..nvars):

  elif test = `TEST_MATRIX`then

    test_array := array(1..nvars,1..nvars):

  fi:



  error_statement := 0:

  error_statement :=  traperror(evalhf(EqProc(xstart, var(test_array))));

  if (type(error_statement,numeric)) then

   RETURN(OK):

  else

lprint(`ERROR in`, convert(EqProc, string)):

lprint(`^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^`);

lprint(``);

lprint(`  ERROR STATEMENT:`,error_statement);

lprint(``);

lprint(`  Error occured in NEWTON/hf/hfNewton when the user defined `);

lprint(`   procedures for the computation of function vector was tested`):

lprint(`  NOTE:`);

lprint(`       A) The head of the procedure must have the form:`):

lprint(`          foo := proc(<STATE_ARRAY>, <RESULT>)`);

lprint(`            ...        `);

lprint(`           end;`);

lprint(`       B) The RESULT is always assigned to the SECOND argument`);

lprint(`          NOTE: if no value is assigned to the second argument or`);

lprint(`                to an element of the array then evalhf returns 0!`);

lprint(`       C) All parameters should be defined GLOBAL. Hence, if you `);

lprint(`          want to use the number of variables in the procedure than`);

lprint(`          define it as a global variable.`);

lprint(`       D) Symbolic packages, like LINALG cannot be used in hardware`);

lprint(`           floating point computations (see MAPLE manual)`);

lprint(``);

   ERROR(error_statement):

  fi:

end:





`NEWTON/hf/ProcTest/ns` := proc(EqProc, xstart, nvars)

  local error_statement,  # error returned by evalhf evaluation

        test_array:       # array for testing the procedure



  test_array := array(1..nvars):

  error_statement := 0:

  error_statement :=  traperror(evalhf(EqProc(xstart, var(test_array),

                                              nvars)));

  if (type(error_statement,numeric)) then

   RETURN(OK):

  else

lprint(`ERROR in`, convert(EqProc, string)):

lprint(`^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^`);

lprint(``);

lprint(`  ERROR STATEMENT:`,error_statement);

lprint(``):

lprint(``);

lprint(`  Error occured in NEWTON/hf/hfNewton when the user defined `);

lprint(`   procedures for the computation of function vector was tested`):

lprint(`  NOTE:`);

lprint(`       A) The head of the procedure must have the form:`):

lprint(`          foo := proc(<UNMODIFIED_DEVIATION>, <STATE_ARRAY>, <NUMBER OF STATES)`);

lprint(`            ...        `);

lprint(`           end;`);

lprint(`       B) The RESULT is always assigned to the SECOND argument`);

lprint(`          NOTE: if no value is assigned to the second argument or`);

lprint(`                to an element of the array then evalhf returns 0!`);

lprint(`       C) All parameters should be defined GLOBAL. Hence, if you `);

lprint(`          want to use the number of variables in the procedure than`);

lprint(`          define it as a global variable.`);

lprint(`       D) Symbolic packages, like LINALG cannot be used in hardware`);

lprint(`           floating point computations (see MAPLE manual)`);

lprint(``);

   ERROR(error_statement):

  fi:

end:



###############################################################

##

##  Procedures for lu-decomposition (taken from BESIRK)

##



`NEWTON/hf/ludecomphf` := proc (A,indx,n)

local i,j,k,tiny,_A,S,Amax,VV,d,dum,imax;

tiny := 1e-20;

d := 1;



VV := array(1..n);

for i from 1 to n do

   Amax := 0;

   for j from 1 to n do

   #####     print(A[i,j]);

     if abs(A[i,j]) > Amax then Amax := abs(A[i,j]) fi;

   od;

   if Amax = 0 then ERROR (`Simgular Matrix`) fi;

   VV[i] := 1 / Amax;

od;

for j from 1 to n do

   for i from 1 to j-1 do

      S := A[i,j];

      for k from 1 to i-1 do

         S := S - A[i,k] * A[k,j]

      od;

      A[i,j] := S;

   od;

   Amax := 0;

   for i from j to n do

      S := A[i,j];

      for k from 1 to j-1 do

         S := S - A[i,k] * A[k,j];

      od;

      A[i,j] := S;

      dum := VV[i] * abs(S);

      if dum >= Amax then imax := i; Amax := dum; fi;

   od;

   if j <> imax then

      for k from 1 to n do

         dum := A[imax,k];

         A[imax,k] := A[j,k];

         A[j,k] := dum;

      od;

      d = - d;

      VV[imax] := VV[j];

   fi;

   indx[j] := imax;

   if A[j,j] = 0 then A[j,j] := tiny; fi;

   if j <> n then

      dum := 1 / A[j,j];

      for i from j+1 to n do

        A[i,j] := A[i,j] * dum;

      od;

   fi;

od;

end:





`NEWTON/hf/lubacksubhf` := proc(A, b, indx, n)

local i,j,k,ii,ll,S;

ii := 0;

for i from 1 to n do

   ll := indx[i];

   S := b[ll];

   b[ll] := b[i];

   if ii <> 0 then

      for j from ii to i-1 do

         S := S - A[i,j] * b[j];

      od;

   elif S <> 0 then

      ii := 1;

   fi;

   b[i] := S;

od;



for i from n to 1 by -1 do

   S := b[i];

   if i < n then

      for j from i+1 to n do

         S := S - A[i,j] * b[j];

      od;

    fi;

   b[i] := S / A[i,i];

od;

end:



###############################################################

##

## Internal procudre computing the jacobian

##   the procedure return the Jacobian matrix in Jmat



`NEWTON/hf/NumericalJacobian` := proc(funcF, ff, xx, JJ, nvars)



local  i, j, k,      # various counters

       ffJ,          # function vector used for the numerical

                     #  computation of the Jacobian

       xxJdel:       # deviation of the state. Used for numerical

                     #  computation of the jacobian

ffJ  := array(1..nvars);



# the Jacobian is evaluated numerically

  for i to nvars do

   # definition of the deviation

   xxJdel := 1e-4*abs(xx[i]);

   if xxJdel < 1e-10 then

    xxJdel := 1e-10:

   fi:

   # definition of the deviation vector and applying it to the

   # function vector

   xx[i] := xx[i]+xxJdel:

   funcF(xx,ffJ):

   # Assignment of the jacobian row:

   for k to nvars do

    if abs(ffJ[k] - ff[k]) < 1e-20 then

     JJ[k,i] := 0

    else

     JJ[k,i] := (ffJ[k] - ff[k])/xxJdel:

    fi:

   od:

   # reset of the deviation vector

   xx[i] := xx[i] - xxJdel:

  # applying procedure to the next variable deviation

  od:

end:



###############################################################

##

## Procedure for iterations in a evalhf environment

##



`NEWTON/hf/NewtonShell` :=

 proc(funcF, funcJ, xstart, nvars, eps, maxtime, maxit, infolevel,

      nextstep, VarPrint)



 local delx,         # change in x vector

       del1,         # adjusted change in x vector

       et,           # Elapsed time

       ff,           # function vector

       F1,           # vector of functions used for numerical Jacobian

       itercount,    # records current iteration number

       i, j, k,      # various counters

       indx,         # [RB03/98] vector for LU decomposition

       JJ,           # jacobian matrix

       nrm,          # norm of function vector

       xnrm,         # norm of variable changes

       starttime,    # Time the calculations started (s)

       timenow,      # Current time (s)

       xstep,        # step size for numerical Jacobian

       xx:           # vector of variable values





# Define a few vectors and matrices for working space

xx   := array(1..nvars):

xx   := xstart;

delx := array(1..nvars);

ff   := array(1..nvars);

F1   := array(1..nvars);

del1 := array(1..nvars);

JJ   := array(1..nvars,1..nvars);





# Evaluate the function vector . The new values are assigned to the

#  arrays ff

funcF(xx,ff);



# Evaluate the jacobian. The new values are assigned to the

#  arrays ff

if not(funcJ = 0) then

 funcJ(xx,JJ);

else

 `NEWTON/hf/NumericalJacobian`(funcF, ff, xx, JJ, nvars);

fi:



#  the norm of the function vector



nrm:=0:

for k to nvars do

 nrm := abs(ff[k]) + nrm:

od:

xnrm := nrm;



# initalization of the permutation vector

indx := array(1..nvars);



# Perform iterations

itercount := 0;

starttime := time();

while (nrm > eps and xnrm > eps) do

   # Print out various things according to inflolevel

   if infolevel[1] > 0 then lprint(`Iteration   `, itercount,

                                  `   Norm = `, nrm);

      print();

   fi;

   if infolevel[2] > 0 then

     lprint(`Variables  `);

     for k to nvars do

      lprint(`    `,`x[`,k,`] = `, xx[k]):

     od:

     print():

   fi;

   if infolevel[2] < 0 then

     lprint(`Variables  `);

     for k to infolevel[7] do

     if VarPrint[k] > 0 then

      lprint(`    `,`x[`,VarPrint[k],`] = `, xx[VarPrint[k]]):

     fi:

     od:

     print():

   fi;

   if infolevel[3] > 0 then

     lprint(`Functions  `);

     for k to nvars do

      lprint(`    `,`f[`,k,`] = `, ff[k]):

     od:

     print():

   fi;

   if infolevel[4] > 0 then

     lprint(`Jacobian  `);

     for k to nvars do

      for j to nvars do

       lprint(`    `,`J[`,k,j,`] = `, JJ[k,j]):

      od

     od:

     print():

   fi;

   # Solve linear system to get vector of corrections

   #     LU decomposition

   #

   #      JJ    get LU decomposed. The function returns a matrix of the

   #             triangular matrices

   #      indx  is the permutation vector

   #      nvars the number of variables

   #

###print-cmd### print(`old x:  `, xx[1], xx[2]):

   evalf(`NEWTON/hf/ludecomphf`(JJ,indx,nvars));

   #     LU re-decomposition and solution of the equation system

   #

   #      JJ    is the LU decomposed matrix

   #      ff    is the left hand side

   #      indx  is the permutation vector of the LU decomposition

   #      nvars the number of variables

   # The result is returned by modifying ff.



   evalf(`NEWTON/hf/lubacksubhf`(JJ,ff,indx,nvars));



   # Compute new set of variable values and the next function vector etc.



   # the nextfunction computes the next state array.



   nextstep(ff,xx,nvars);



  # Evaluate the function vector and jacobian matrix. The new

  #  values are assigned to the arrays ff and JJ



  funcF(xx,ff);



  nrm:=0:

  for k to nvars do

   nrm := abs(ff[k]) + nrm:

  od:



  # Evaluate the jacobian. The new values are assigned to the

  #  arrays ff

  if not(funcJ = 0) then

   funcJ(xx,JJ);

  else

   `NEWTON/hf/NumericalJacobian`(funcF, ff, xx, JJ, nvars);

  fi:



  #  the norm of the variable change vector



  xnrm := 0;

  for k to nvars do

   xnrm := abs(ff[k]) + xnrm:

  od:



  # Increment iteration count and check that we have not done too

  #  many iterations

  itercount := itercount+1;

  if itercount > (maxit + 1) then

   ERROR(`Maximum number of iterations exceeded`);  fi;



  # Check time

  timenow := time();

  et := timenow - starttime;

  if et > maxtime then ERROR(`Time limit exceeded`); fi;

od:



# Returning the result

for k to nvars do

 xstart[k] := xx[k];

od:

end: