2D Newton's method in python

Question

I am working on coding a backward Euler method in Python and I am having problems coding the Newton part. We are given a tolerance of 1e-4 and using this I am getting very small numbers in the output vector for my Newton's method.

This is my code:

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)

For this problem, I am using these function functions:

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)

When I run it I get a matrix:

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

which is way too small to use in my backwards Euler function, which makes me think there is something wrong with my Newton. Can anyone help me figure out what I am doing wrong here? From my understanding of the Newton function this seems like it should be the correct thing, but clearly it's not doing exactly what I need it to do.

Also to run this function at the top of my code I have:

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

which are not all needed for this, but some are!

I am just trying to get the best newton approximation so I can use my backwards euler function. The problem with the Newton output I get is that I would need to do thousands of BE iterations in order to get any answer, which shouldn't happen given the accuracy ability of the BE function. — bananagurlz, Jun 11 '18 at 20:54
I can also add my Backward Euler code if it helps make this problem make more sense — bananagurlz, Jun 11 '18 at 21:01
By the way, `error=xn-x=(x-y0)-x=-y0`. I would rename `y0` as `correction` and check that correction is smaller than `tol` in the `while` loop. — AGN Gazer, Jun 11 '18 at 21:15

score 0 · Answer 1 · answered Jun 11 '18 at 21:07

0

There may not be a problem with your solution (except for the computation of the Jacobian). It seems to me that you are trying to solve A*x==0. The obvious solution to this system is x=0 which is what you get.

To see this, try changing your function() to: z = A * x + np.matrix([[1], [-3]])

answered Jun 11 '18 at 21:07

AGN Gazer

8,025
2
27
45

Can you explain why I would do this one more time? Sorry I am still confused. Where does the second matrix come from? Is this a problem in the functional computation of the Jacobian or the actual function itself? if it is the function, that specific code for the function worked for other iterative methods, so what is different now? – bananagurlz Jun 11 '18 at 21:26
I tried doing this and it gave me much nicer answers, thank you! But I am still unsure why. Also now when I set my newton value as xn and then run my new function (A,xn) it gives me the zero vector. Is this for the same reason? Or are we solving inverses somehow? – bananagurlz Jun 11 '18 at 21:35
@bananagurlz All I did was to change your function (or system of equations) so that the solution if not a trivial one. Let's think in terms of 1D case. What _you_are doing is equivalent to solving `-33.4*x==0`. What is the solution of this equation? Obviously, it is 0. What did I do? I changed the equation to something like `-33.4*x +1 == 0`. Now the solution is different from zero. So, the issue (I think) is that your system of equations (as defined by you) has trivial `(0, 0)` solution. You need to change your equations if you want to find non-trivial solutions. – AGN Gazer Jun 11 '18 at 22:07
That makes sense. So then I need find a way to write a program to solve this system. Was the 1,3 vector arbitrary or did you have a reason for it? If you did, how did you solve it? I am just trying to think about how it would be solved. – bananagurlz Jun 11 '18 at 22:11
@bananagurlz It was arbitrary just to illustrate the point. You need to re-think your problem if you have to. I cannot tell you what system you need to solve. – AGN Gazer Jun 11 '18 at 22:13
Ok. Just for a hint, would I be using the matrix inverse somehow? I guess I am still confused why A*x where x is clearly a vector is not the right way to solve it. Would I have A*x=y (y unknown) and do like $A^{-1} y=x$ so then it solves an actual system? – bananagurlz Jun 11 '18 at 22:18
@bananagurlz A simple linear system of equations such as `A*x==y` can be solved using matrix inverse directly. You do not need Newton method for this. So, the thing is that it is _you_ who needs to know what you are solving. Is it A*x==0? Is it `f1(x,y) == 0; f2(x, y)==0` where `f1` and `f2` are non-lienear functions of `x` and `y`? I cannot tell you what your problem is that you are trying to solve. – AGN Gazer Jun 11 '18 at 22:26
@bananagurlz Take a look at https://math.stackexchange.com/q/1544748 and other helpful resources online. – AGN Gazer Jun 11 '18 at 22:28
I am sorry if I am not making a ton of sense. What my newton is, given x0 and A, it is finding a solution for xn. It is using the formula x(n+1)=xn-(f(xn)/f'(xn)). Basically it should solve out enough xn's to get an approximation close enough to the actuall solution to the equation, f(x). So basically I am looking for it to find an xn where it is within 1e-4 distance between itself and the previous xn. This might have been something you already knew. I am not trying to make you talk in circles and I hope I am not disturbing you, but I just want to put myself on the right track. – bananagurlz Jun 11 '18 at 22:32
@bananagurlz `x0` and `xn` are basically the same thing: first one is initial approximation to the solution while `xn` is the n-th approximation to the solution of the same system of equations. What I am saying is that you are solving the equation `function=A*x == 0`. I can tell you right away that if `det(A)!=0` the solution is `x=[0,0]`. You do not need Newton method at all - the solution is bvious. – AGN Gazer Jun 11 '18 at 22:36
@bananagurlz I think you did not understood the problem well. I think you are trying to solve equation (10) in http://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node3.html – AGN Gazer Jun 11 '18 at 22:53
that is what I am trying to solve, I understand it in one dimension, and now understand from what you have said why it is giving me an error (thank you!), it is just going to be a matter of finding out how to fix it. – bananagurlz Jun 11 '18 at 23:10

2D Newton's method in python

1 Answers1