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);
   7     FuncList  - list of functions to solve.
   8     VarList   - list of variable names.
   9     GuessList - list of initial approximations.
  10
  11 Flags:
  12   newtonepsilon
  13     default: [10.0^(-fpprec/2)] - precision to determine when the
  14     process has converged towards the solution.
  15   newtonmaxiter
  16     default: [50] - maximum number of iterations to stop the process
  17     if this one does not converge or if it converges too much slowly.
  18
  19 Results:
  20   The solution is returned in the same format that SOLVE() returns.
  21   If the solution isn't found, [] is returned.
  22
  23 History:
  24   2003-11 Salvador Bosch Perez - version 1.0
  25
  26 --
  27
  28 Copyright (C) 2003  Salvador Bosch Perez
  29
  30 This library is free software; you can redistribute it and/or
  31 modify it under the terms of the GNU Lesser General Public
  32 License as published by the Free Software Foundation; either
  33 version 2.1 of the License, or (at your option) any later version.
  34
  35 This library is distributed in the hope that it will be useful,
  36 but WITHOUT ANY WARRANTY; without even the implied warranty of
  37 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  38 Lesser General Public License for more details.
  39
  40 You should have received a copy of the GNU Lesser General Public
  41 License along with this library; if not, write to the Free Software
  42 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  43 */
  44
  45 newtonepsilon:10.0^(-16/2)$ /* Default vaule is 10.0^-8 */
  46 newtonmaxiter:50$
  47 newtondebug:false$
  48
  49 solve_by_lu(eqs, vars, field) := block(
  50     [M, b, n : length(vars), sol],
  51     M : augcoefmatrix(eqs, vars),
  52     b : -col(M, n+1),
  53     M : submatrix(M, n+1),
  54     sol : linsolve_by_lu(M, b, field),
  55     map("=", vars, first(args(transpose(sol[1]))))
  56 )$
  57
  58 mnewton(FuncList, VarList, GuessList):=
  59   block([nfunc, Solutions, Increments, solved:false, h, DF, numer:numer,
  60          i, j, k, keepfloat:true, ratprint:false, field : 'complexfield, float2bf : true],
  61     local(h),
  62
  63     if (not(listp(FuncList))) then FuncList: [FuncList],
  64     if (not(listp(VarList))) then VarList: [VarList],
  65     if (not(listp(GuessList))) then GuessList: [GuessList],
  66
  67     /* Decide which field to use */
  68     if bfloatp(newtonepsilon) then field : 'bigfloatfield,
  69
  70     if field = 'bigfloatfield then
  71        FuncList : map(lambda([u], lhs(u)-rhs(u)), bfloat(FuncList))
  72     else (
  73        numer : true,
  74        FuncList : map(lambda([u], lhs(u)-rhs(u)), float(FuncList))),
  75
  76     /* Depuration of the function arguments */
  77     nfunc:length(FuncList),
  78     if length(VarList) # nfunc then (
  79       print("mnewton: incorrect number of variable names (", nfunc,
  80             "functions but", length(VarList), "variable names)."),
  81       return(false)
  82     ),
  83     if length(GuessList) # nfunc then (
  84       print("mnewton: incorrect number of approach values (", nfunc,
  85             "variables but", length(GuessList), "approximation values)."),
  86       return(false)
  87     ),
  88     /* apply(kill, VarList), */
  89
  90     DF: zeromatrix (nfunc, nfunc),
  91     h : makelist (h[i], i, 1, nfunc),
  92     for i:1 thru nfunc do
  93       for j:1 thru nfunc do
  94         DF[i][j] : diff (FuncList[i], VarList[j]),  /* Jacobian matrix */
  95
  96     if field = 'complexfield then
  97       GuessList: float (GuessList)
  98     else
  99       GuessList: bfloat(GuessList),
 100
 101     if newtondebug
 102        then print ("DEBUG: initial GuessList:", GuessList),
 103
 104     /* Solution of the functions system */
 105     for k:1 thru newtonmaxiter do (
 106       Solutions:map("=", VarList, GuessList),
 107
 108       /* Solve 0 = F(x) + DF(x) h for h.
 109        */
 110       block
 111         ([DFx: subst (Solutions, DF),
 112         Fx: subst (Solutions, FuncList)],
 113         Increments: solve_by_lu (transpose (DFx . h + Fx)[1], h, field)),
 114
 115       GuessList:GuessList + map ('rhs, Increments),
 116       if field = 'complexfield then
 117         GuessList: float (GuessList)
 118       else
 119         GuessList: bfloat(GuessList),
 120
 121       if newtondebug
 122          then print ("DEBUG: Increments:", Increments, "GuessList:", GuessList),
 123
 124       solved:true,
 125       for i:1 thru nfunc do
 126         solved: solved and abs (rhs (Increments[i])) < newtonepsilon,
 127       if solved then return(0)
 128     ),
 129
 130     /* Return of solution or error */
 131     if solved = false then (
 132       print("mnewton: the process doesn't converge or it converges too slowly."),
 133       return([])
 134     ),
 135     Solutions:map("=", VarList, GuessList),
 136     if newtondebug
 137        then print ("DEBUG: subst (Solutions, FuncList) =>", subst (Solutions, FuncList)),
 138     return([Solutions])
 139   )$