share/mnewton/mnewton.mac

   1 /*  mnewton.mac v1.1 for Maxima (tested with Maxima 5.9.0).
   2
   3 Multiple nonlinear functions solution using the Newton method.
   4
   5 Usage:
   6   mnewton(FuncList, VarList, GuessList, [DF]);
   7     FuncList  - list of functions to solve.
   8     VarList   - list of variable names.
   9     GuessList - list of initial approximations.
  10     DF        - Jacobian matrix (optional)
  11
  12 Flags:
  13   newtonepsilon
  14     default: [10.0^(-fpprec/2)] - precision to determine when the
  15     process has converged towards the solution.
  16   newtonmaxiter
  17     default: [50] - maximum number of iterations to stop the process
  18     if this one does not converge or if it converges too much slowly.
  19
  20 Results:
  21   The solution is returned in the same format that SOLVE() returns.
  22   If the solution isn't found, [] is returned.
  23
  24 History:
  25   2003-11 Salvador Bosch Perez - version 1.0
  26
  27 --
  28
  29 Copyright (C) 2003  Salvador Bosch Perez
  30
  31 This library is free software; you can redistribute it and/or
  32 modify it under the terms of the GNU Lesser General Public
  33 License as published by the Free Software Foundation; either
  34 version 2.1 of the License, or (at your option) any later version.
  35
  36 This library is distributed in the hope that it will be useful,
  37 but WITHOUT ANY WARRANTY; without even the implied warranty of
  38 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  39 Lesser General Public License for more details.
  40
  41 You should have received a copy of the GNU Lesser General Public
  42 License along with this library; if not, write to the Free Software
  43 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  44 */
  45
  46 newtonepsilon:10.0^(-16/2)$ /* Default value is 10.0^-8 */
  47 newtonmaxiter:50$
  48 newtondebug:false$
  49
  50 mnewton(FuncList, VarList, GuessList, [DF]):=
  51   block([nfunc, Solutions, Increments, solved:false, extra_vars,
  52          i, j, k, ratprint:false, field : 'complexfield, float2bf : true],
  53
  54     if (not(listp(FuncList))) then FuncList: [FuncList],
  55     if (not(listp(VarList))) then VarList: [VarList],
  56     if (not(listp(GuessList))) then GuessList: [GuessList],
  57
  58     /* Refuse to proceed if FuncList depends on some variable other than VarList.
  59      * Even in that case (which means that the result will just be an expression
  60      * containing any other variables), it is plausible that mnewton could carry
  61      * out a fixed number of iterations, but mnewton would have to reorganized
  62      * somewhat to implement that. For now require that FuncList depends only on VarList.
  63      */
  64
  65     if (extra_vars: setdifference (setify (listofvars (FuncList)), setify (VarList))) # {}
  66       then (printf (true, "mnewton: extra variables in equations; found: ~{~a~^, ~}~%", extra_vars),
  67             return (false)),
  68
  69     /* Decide which field to use */
  70     if bfloatp(newtonepsilon) then field : 'bigfloatfield,
  71
  72     if field = 'bigfloatfield then
  73        FuncList : map(lambda([u], lhs(u)-rhs(u)), bfloat(FuncList))
  74     else (
  75        FuncList : map(lambda([u], lhs(u)-rhs(u)), float(FuncList))),
  76
  77     /* Depuration of the function arguments */
  78     nfunc:length(FuncList),
  79     if length(VarList) # nfunc then (
  80       print("mnewton: incorrect number of variable names (", nfunc,
  81             "functions but", length(VarList), "variable names)."),
  82       return(false)
  83     ),
  84     if length(GuessList) # nfunc then (
  85       print("mnewton: incorrect number of approach values (", nfunc,
  86             "variables but", length(GuessList), "approximation values)."),
  87       return(false)
  88     ),
  89     /* apply(kill, VarList), */
  90
  91     /* Jacobian matrix */
  92     if DF = []
  93       then (DF: zeromatrix (nfunc, nfunc),
  94             for i:1 thru nfunc do
  95               for j:1 thru nfunc do
  96                 DF[i][j] : diff (FuncList[i], VarList[j]))
  97        else DF: DF[1],
  98
  99     if newtondebug then print ("DEBUG: DF =", DF),
 100
 101     if field = 'complexfield then
 102       GuessList: float (GuessList)
 103     else
 104       GuessList: bfloat(GuessList),
 105
 106     if newtondebug
 107        then print ("DEBUG: initial GuessList:", GuessList),
 108
 109     /* Solution of the functions system */
 110     for k:1 thru newtonmaxiter do (
 111       Solutions:map("=", VarList, GuessList),
 112
 113       /* Solve 0 = F(x) + DF(x) h for the increment h.
 114        */
 115       block
 116         ([DFx: subst (Solutions, DF),
 117         Fx: subst (Solutions, FuncList),
 118         linsolve_solution],
 119         linsolve_solution: linsolve_by_lu (DFx, - Fx, field),
 120         Increments: first (args (transpose (linsolve_solution[1])))),
 121
 122       GuessList: GuessList + Increments,
 123       if field = 'complexfield then
 124         GuessList: float (GuessList)
 125       else
 126         GuessList: bfloat(GuessList),
 127
 128       if newtondebug
 129          then print ("DEBUG: Increments:", Increments, "GuessList:", GuessList),
 130
 131       solved:true,
 132       for i:1 thru nfunc do
 133         solved: solved and abs (Increments[i]) < newtonepsilon,
 134       if solved then return(0)
 135     ),
 136
 137     /* Return of solution or error */
 138     if solved = false then (
 139       print("mnewton: the process doesn't converge or it converges too slowly."),
 140       return([])
 141     ),
 142     Solutions:map("=", VarList, GuessList),
 143     if newtondebug
 144        then (print ("DEBUG: Solutions =", Solutions),
 145              print ("DEBUG: FuncList =", FuncList),
 146              print ("DEBUG: subst (Solutions, FuncList) =>", subst (Solutions, FuncList))),
 147
 148     return([Solutions])
 149   )$