我正试图创建一个函数,为多个变量(dx=-Jinv^-1*F(进行Newton-Raphson寻根,但它不会连续迭代一个循环(即输出不会成为下一次迭代的下一个输入(。有人看到这个问题了吗?
PS。我试图在不使用任何内置矩阵函数的情况下对其进行编码,因此这些解决方案将不适用于此。
# Multivariable Root Finding via Newton Raphson
# import necessary functions
import numpy as np
# define array function for f1 and f2
x = np.arange(-3,3,.01)
def f1(x):
return 2*(x[0]**3) - 6*(x[0]*(x[1]**2)) - 1
def f2(x):
return -2*(x[1]**3) + 6*(x[1]*(x[0]**2)) + 3
# define matrix F(x) 2x1
def F(x_vec):
return [f1(x_vec),f2(x_vec)]
#print(F([2,3]))
# Create Inverse Jacobian 2x2
# InvJacobian Upper Left
def JUL(x):
return (6*x[0]**2 - 6*x[1]**2)/((6*x[0]**2 - 6*x[1]**2)**2 + 144*(x[0]**2)*(x[1]**2))
# InvJacobian Upper Right
def JUR(x):
return (12*x[0]*x[1])/((6*x[0]**2 - 6*x[1]**2)**2 + 144*(x[0]**2)*(x[1]**2))
# InvJacobian Lower Left
def JLL(x):
return (-12*x[0]*x[1])/((6*x[0]**2 - 6*x[1]**2)**2 + 144*(x[0]**2)*(x[1]**2))
# InvJacobian Lower Left
def JLR(x):
return (6*x[0]**2 - 6*x[1]**2)/((6*x[0]**2 - 6*x[1]**2)**2 + 144*(x[0]**2)*(x[1]**2))
# Combine all Jacobian into Matrix
def Jinv(x_vec):
return [JUL(x_vec),JUR(x_vec),JLL(x_vec),JLR(x_vec)]
# Create Newton Raphson Functon
def rf_newton2d(F_system, Jinv_system, x_vec0, tol, maxiter):
# let ou be the starting array x_vec0
ou = x_vec0
i = 0
err = 10**-4
# F*Jinv should output a 2x1 matrix
while err > tol and i <= maxiter:
F_system = F(ou)
Jinv_system = Jinv(ou)
# w is the [0] position of the 2x1 output matrix
w = F_system[0]*Jinv_system[0] + F_system[1]*Jinv_system[1]
# z is the [1] position of the 2x1 output
z = F_system[0]*Jinv_system[2] + F_system[1]*Jinv_system[3]
# u is the output matrix of F*Jinv
u = [w,z]
# out is the new array involving u + ou
out = [ou[0] + u[0],ou[1] +u[1]]
#define nu as out
nu = out
# let nu be the ou for the next iteration
ou = nu
#run continuous iterations
i += 1
return nu
print(rf_newton2d(F([x_vec0]),Jinv([x_vec0]),[2,3],10**-5, 5))
return
语句会立即停止该函数,并从调用该函数的位置继续运行主代码。在我看来,while
循环中的return
语句不应该在while
循环中。
在最后一行中,您有函数F([x_vec0])
,这意味着函数F
中有变量x_vec = [x_vec0]
。然后将此变量传递给函数f1
,因此在f1
中为x = [x_vec0]
。在x[1]
内可以在f1
内找到。问题是x
永远不会有第二个条目。如果x_vec0 = [0,1]
(或一些其他数组对象,如np.ndarray
(,则不存在[x_vec0] = [[0,1]]
,即[x_vec0][0] = [0,1]
和[x_vec0][1]
。
我从回溯中得到了所有信息(显然我必须为x_vec0
创造一个值(:
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-4-d57a9d10dd94> in <module>
80 return nu
81 x_vec0 = np.array([1,2])
---> 82 print(rf_newton2d(F([x_vec0]),Jinv([x_vec0]),[2,3],10**-5, 5))
<ipython-input-4-d57a9d10dd94> in F(x_vec)
16 # define matrix F(x) 2x1
17 def F(x_vec):
---> 18 return [f1(x_vec),f2(x_vec)]
19
20 #print(F([2,3]))
<ipython-input-4-d57a9d10dd94> in f1(x)
9
10 def f1(x):
---> 11 return 2*(x[0]**3) - 6*(x[0]*(x[1]**2)) - 1
12
13 def f2(x):
IndexError: list index out of range
如果不知道牛顿-拉斐森根的数学原理,我就无法提出解决方案。