python中的二维牛顿法

Question

我正在编写 Python 中的后向 Euler 方法，但在编写 Newton 部分时遇到问题。我们得到了 1e-4 的容差，使用这个我在牛顿法的输出向量中得到了非常小的数字。

这是我的代码：

def Newt(U,x,tol): #takes in matrix, x vector, tolerance 

    error=np.matrix([[1],[1]])  #set an error matrix of 1,1
    while abs(LA.norm(error))>tol:
        func=function(U,x)   #define a f(x) vector
        jac=functionprime(U,x) #define the inverse jacobian vector
        y0=jac*func  #the change vector is the jacobian inverse times the function
        xn=x-y0 #the new x is the difference
        error=xn-x #set the error
        x=xn
    print(x)     

对于这个问题，我使用了这些函数函数：

A=np.matrix([[-33.4, 66.6], [33.3, -66.7]]) #A matrix
x0=np.matrix([[3],[0]]) #x0 verctor

def function(A,x):
    z=A*x
    return(z) #all my function does is multiply the matrix by the x vector

def functionprime(A,x0):
    b=(1/10)*np.matrix([[-66.7,-66.6],[-33.3,-33.4]]) #tried to code this to just output the inverse jacobian 
    return(b)

当我运行它时，我得到一个矩阵：

[[  4.31664687e-27],
 [  2.15832344e-27]]

这太小了，无法在我的反向欧拉函数中使用，这让我觉得我的牛顿有问题。谁能帮我弄清楚我在这里做错了什么？根据我对牛顿函数的理解，这似乎应该是正确的，但显然它并没有完全按照我的需要做。

还要在我的代码顶部运行这个函数：

import matplotlib.pyplot as plt
import math
import numpy as np
from pylab import *
from numpy import linalg as LA

这不是全部都需要，但有些是！

Answer 1

您的解决方案可能没有问题（雅可比计算除外）。在我看来，您正在尝试解决A*x==0。这个系统的明显解决方案是x=0，这就是你得到的。

要查看此内容，请尝试将您的 function() 更改为：z = A * x + np.matrix([[1], [-3]])