I've written a code in python which implements the Newton-Raphson method to solve multiple nonlinear equations.
The specific question I've taken is from Mark Newman's - Computational Physics, exercise 6.17 Nonlinear circuits
import numpy as np
from numpy.linalg import solve, norm
from math import e
#DATA
vp= 5. #V_plus in volts
r1, r2, r3, r4 = 1., 4., 3., 2. #k-ohm resistances
i = 3. #A constant originally in nA
vt = 0.05 #V_t in volts
def f(x):
'''
evaluates f(x) for where x is a 2-dim vector of voltage v1 and v2
'''
v1, v2 = x[0], x[1]
y = np.array([(v1-vp)/r1 + v1/r2 + i*(e**((v1-v2)/vt)-1), (vp-v2)/r3 - v2/r4 + i*(e**((v1-v2)/vt)-1)])
print y
return y
def gradf(x):
'''
n-Derivative of f(x) where x is a vector of n-dimensions
'''
v1, v2 = x[0], x[1]
m = np.array([[1./r1 + 1./r2 + (i/vt)*e**((v1-v2)/vt), (i/vt)*e**((v1-v2)/vt)],\
[(-i/vt)*e**((v1-v2)/vt), -1*(1./r3 +1./r4 +(i/vt)*e**((v1-v2)/vt))]], dtype = np.float64)#the matrix for the 'grad' f function
print m
return m
def cls_newton(x):
'''
Classroom implementation of the newton raphson method
'''
v1, v2 = x[0], x[1]
f_v1 = 1./r1 + 1./r2 + (i/vt)*e**((v1-v2)/vt)
f_v2 = (-i/vt)*e**((v1-v2)/vt)
g_v1 = (i/vt)*e**((v1-v2)/vt)
g_v2 = -1*(1./r3 +1./r4 +(i/vt)*e**((v1-v2)/vt))
f = (v1-vp)/r1 + v1/r2 + i*(e**((v1-v2)/vt)-1)
g = (vp-v2)/r3 - v2/r4 + i*(e**((v1-v2)/vt)-1)
print f
print g
print f_v2, g_v1, g_v1, f_v1
v1n = v1 - (f*g_v2 - g*f_v2)/(f_v1*g_v2 - g_v1*f_v2)
v2n = v2 - (g*f_v1 - f*g_v1)/(f_v1*g_v2 - g_v1*f_v2)
print v1n
print v2n
return np.array([v1n,v2n])
x1 = np.array([4., 5.]) #initial guess of roots are 4. and 5. volts
error = 1e-6 # permissible error
i = 0 # iteration counter
while norm(x1)>error and i < 50:
delta = solve(gradf(x1), f(x1))
x2 = x1 - delta
print x2
print 'x1 = {0}, x2 = {1}'.format(x1, x2)# test line
x1 = x2
print x1
i+=1
rt = x1 # estimated root of the equation
print 'The root of the equation is ' + str(rt) + '\n' + 'f(root) = ' + str(f(rt))
print 'No. of iterations: ' + str(i)
In this code I've written functions for two different implementations of the method for multiple roots.
The one which I've used in this program is one where I solve an equation between gradf(x) (Which produces a Jacobian matrix) and f(x) (Which gives me the a vector with the equations which I found using Kirchoff's laws).
it works like gradf(x).delta = f(x) so we find delta using the solve() function and we subtract delta from x1 (our initial v1 and v2) to find x2
I'm having a problem with the matrix though, when I call the function gradf([4.,5.]) in Ipython, it gives me a matrix like
array([[ 1.25000004e+00, 4.12230724e-08],
[ -4.12230724e-08, -8.33333375e-01]])
but the same matrix when printed during the normal operation of the program is something like
[[ 1.25 0. ]
[ 0. -0.83333333]]
I get this matrix in my first iteration regardless of the initial guess of v1 and v2 (or x1). The next iteration gives me an error like
LinAlgError: Singular matrix .
I don't think this is due to rounding off in Python either because when I individually print the value of (say) the first array element in the matrix (while running the script), it gives me a zero where it should be giving something like 4.12230724e-08.
The classroom implementation or cls_newton(x) which simplifies the equations before-hand and directly gives me x2 seems to do the same thing but I can't tell why, it gives me a different answer through Ipython and a different one during execution.
Also, when I write say f_v1 I'm referring to the partial derivative of f with respect to v1 and g_v2 would be the partial derivative of g with respect to v2 and so on.
Thank you for your help in advance!
LinAlgError? A single element won't tell you much about the determinant, which tells you if the matrix is singular. - Andras Deak