2D метод Ньютона в питоне - PullRequest
Новая
       7

2D метод Ньютона в питоне

0 голосов
/ 11 июня 2018

Я работаю над кодированием обратного метода Эйлера в Python, и у меня возникают проблемы с кодированием части Ньютона.Нам дают допуск 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

, которые не все необходимы для этого, но некоторые из них!

1 Ответ

0 голосов
/ 12 июня 2018

Может не быть проблемы с вашим решением (за исключением вычисления якобиана).Мне кажется, что вы пытаетесь решить A*x==0.Очевидным решением для этой системы является x=0, что вы и получаете.

Чтобы увидеть это, попробуйте изменить function() на: z = A * x + np.matrix([[1], [-3]])

...