Attribute VB_Name = "NROOTS"
' Copyright © 2019 John S. Philo
' This file is part of the program SEDNTERP 1. SEDNTERP 1 is free software: you can redistribute it and/or modify
' it under the terms of the GNU General Public License as published by the Free Software Foundation,
' either version 3 Of the License, or any later version. The author asks only that if you use any of
' this code in your own programs then please acknowledge that use and cite any appropriate publications.
' This program is distributed in the hope that it will be useful,
' but WITHOUT ANY WARRANTY; without even the implied warranty of
' MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
' GNU General Public License for more details.
' You should have received a copy of the GNU General Public License
' along with this program (the file Copying.txt). If not, see
'
' these routines were converted to Basic from Press's Numerical Methods
' these routines require
' Sub UsrFun (guess() As Double, alpha1() As Double, beta() As Double)
' where guess is the current solution vector
' alpha1 is a matrix of partial derivatives,
' beta is the negative of the function values
DefInt I-P
DefDbl A-H, R-Z
Sub LUBacksub(A() As Double, N As Integer, np As Integer, indx() As Integer, B() As Double)
' A() is input n x n matrix from LUDecomp, with dimension np x np
' B() is input right-hand side vector, dimension n
' B() is output with solution vector
ii = 0
' when ii becomes positive it will point to first non-zero element of B()
For I = 1 To N
ll = indx(I)
Sum = B(ll)
B(ll) = B(I)
If ii <> 0 Then
For j = ii To I - 1
Sum = Sum - A(I, j) * B(j)
Next j
ElseIf Sum <> 0 Then
ii = I
End If
B(I) = Sum
Next I
' now do the backsubstitution
For I = N To 1 Step -1
Sum = B(I)
If I < N Then
For j = I + 1 To N
Sum = Sum - A(I, j) * B(j)
Next j
End If
B(I) = Sum / A(I, I)
Next I
End Sub
Sub LUDecomp(A() As Double, N As Integer, np As Integer, indx() As Integer, d As Integer)
' A() is input n x n matrix, with dimension np x np
' A() is output with rearrangements recorded in indx()
' d is output as +/- 1 to indicate even or odd rearrangements
Const tiny = 1E-37
d = 1
ReDim VV(1 To np)
For I = 1 To N
aamax = 0
For j = 1 To N
If Abs(A(I, j)) > aamax Then aamax = Abs(A(I, j))
Next j
If aamax = 0 Then
MsgBox "Singular matrix, aborting...", 32, TheProgram
ESCFlag = True
Exit Sub
End If
VV(I) = 1# / aamax ' save scaling
Next I
For j = 1 To N ' loop over columns of Crout's method
If j > 1 Then
For I = 1 To j - 1
Sum = A(I, j)
If I > 1 Then
For k = 1 To I - 1
Sum = Sum - A(I, k) * A(k, j)
Next k
A(I, j) = Sum
End If
Next I
End If
aamax = 0
For I = j To N
Sum = A(I, j)
If j > 1 Then
For k = 1 To j - 1
Sum = Sum - A(I, k) * A(k, j)
Next k
A(I, j) = Sum
End If
dum = VV(I) * Abs(Sum) ' figure of merit for pivot
If dum > aamax Then ' check for best so far
imax = I
aamax = dum
End If
Next I
If j <> imax Then ' we need to interchange rows
For k = 1 To N
dum = A(imax, k)
A(imax, k) = A(j, k)
A(j, k) = dum
Next k
d = -d ' change parity of d
VV(imax) = VV(j) ' also interchange scale factor
End If
indx(j) = imax
If j <> N Then ' divide by the pivot element
If A(j, j) = 0# Then A(j, j) = tiny ' it may be better to be tiny than zero
dum = 1# / A(j, j)
For I = j + 1 To N
A(I, j) = A(I, j) * dum
Next I
End If
Next j ' back to next column
If A(N, N) = 0# Then A(N, N) = tiny
End Sub
Sub MNewton(ntrial As Integer, guess() As Double, N As Integer, tolx As Double, tolf As Double)
' ntrial is max steps
' guess() contains initial guesses on entry, refined on exit
' n is dimension of guess
' tolx is tolerance limit for sum of absolute increments of variables
' tolf is summed absolute function deviations from zero
' assumes routine UsrFun will return function and derivatives
ReDim indx(1 To N), alpha1(1 To N, 1 To N), beta(1 To N)
For k% = 1 To ntrial
UsrFun guess(), alpha1(), beta()
errf = 0
' check function convergence
For I = 1 To N
errf = errf + Abs(beta(I))
Next I
If errf < tolf Then Exit For
' solve using LU decomposition
LUDecomp alpha1(), N, N, indx(), d%
LUBacksub alpha1(), N, N, indx(), beta()
' check root convergence and update solution
errx = 0
For I = 1 To N
errx = errx + Abs(beta(I))
guess(I) = guess(I) + beta(I)
Next I
If errx < tolx Then Exit For
Next k%
If k% = ntrial Then
' bad news
End If
End Sub