牛顿-拉夫森法python实现中的奇异矩阵

Question

我用python编写了一个代码，它实现了求解多个非线性方程的牛顿-拉夫森方法。

我回答的具体问题来自马克·纽曼的“计算物理”，练习6.17非线性电路

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)

在这段代码中，我为多根方法的两个不同实现编写了函数。

我在这个程序中使用的方法是求解梯度f(X)(产生Jacobian矩阵)和f(x)之间的方程(它给我一个向量，其中包含我用Kirchoff定律找到的方程)。

它的工作方式类似于gradf(x).delta = f(x)，因此我们使用so ()函数找到增量，并从x1 (我们的初始v1和v2)中减去增量以找到x2。

但是，当我调用函数梯度f(4，5)时，我对矩阵有一个问题。在Ipython中，它给出了一个矩阵，如

array([[  1.25000004e+00,   4.12230724e-08],
   [ -4.12230724e-08,  -8.33333375e-01]])

但是，在程序正常运行期间打印相同的矩阵时，如下所示

[[ 1.25        0.        ]
[ 0.         -0.83333333]]

我在第一次迭代中得到了这个矩阵，而不考虑v1和v2 (或x1)的最初猜测。下一次迭代给出了一个错误，如

LinAlgError: Singular matrix .

我不认为这是由于在Python中舍入，因为当我单独打印矩阵中的第一个数组元素的值(比如说)时(运行脚本)，它给出了一个零，它应该给出类似于4.12230724e-08的值。

--课堂实现或cls_newton(x) --它简化了前面的公式，直接给了x2 --似乎也在做同样的事情，但我不知道为什么，它通过Ipython给了我一个不同的答案，在执行过程中给出了不同的答案。

，在我写f_v1时，我指的是f相对于v1的偏导数，g_v2是g相对于v2的偏导数，等等。

谢谢你提前提供帮助！

Answer 1

我知道这是一个很晚的反应。如果您仍然对得到答案感兴趣，我认为这段代码应该能做到这一点。对于问题的最后一部分，它给了我正确的答案。如果你有任何问题，请告诉我。

import numpy as np

Vp = 5 #V
R1 = 1e3 # ohms
R2 = 4e3 # ohms
R3 = 3e3 # ohms
R4 = 2e3 # ohms
I0 = 3e-9 # A
VT = 0.05 # V

tol = 1e-6

def f1(V):
    v1 = V[0]; v2 = V[1]
    return( (v1-Vp)/R1 + v1/R2 + I0*(np.exp((v1-v2)/VT)-1) )

def f2(V):
    v1 = V[0]; v2 = V[1]
    return( -(v2-Vp)/R3 - v2/R4 + I0*(np.exp((v1-v2)/VT)-1) )

def j11(V):
    v1 = V[0]; v2 = V[1]
    return 1/R1 + 1/R2 + I0/VT * np.exp((v1-v2)/VT)

def j12(V):
    v1 = V[0]; v2 = V[1]
    return -I0/VT * np.exp((v1-v2)/VT)

def j21(V):
    v1 = V[0]; v2 = V[1]
    return I0/VT * np.exp((v1-v2)/VT)

def j22(V):
    v1 = V[0]; v2 = V[1]
    return -1/R3 - 1/R4 - I0/VT * np.exp((v1-v2)/VT)

# initial guesses
v1 = 0.5
v2 = 0.5

V = np.array( [v1,v2] )
F = np.array( [ f1( V ) , f2( V ) ] )
J = np.array( [ [ j11( V ) , j12( V ) ] , [ j21( V ) , j22( V )] ] )
DV = np.dot( np.linalg.inv(J) , F )
estimate = V - DV
err = np.abs(estimate - V)
while ( err > tol ).any():
    F = np.array( [ f1(estimate) , f2(estimate) ] )
    J = np.array( [ [ j11(estimate) , j12(estimate) ] , [ j21(estimate) , j22(estimate) ] ] )
    DV = np.dot( np.linalg.inv(J) , F ) # f(x)/f'(x)
    new_estimate = estimate - DV
    err = np.abs(new_estimate - estimate)
    estimate = new_estimate

print("V1={:.4E}V\tV2={:.4E}V".format(*new_estimate))
print("Voltage across forward biased diode: {:.4E}V".format(new_estimate[0]-new_estimate[1]))