[Python] Newton-Raphson法 収束するまでの途中経過の描画

import random
import numpy as np
import math
import matplotlib.pyplot as plt

def f1(x1,x2):
    return x2**3-3*x1**2*x2

def f2(x1,x2):
    return x1**3-3*x1*x2**2-4

def df1_dx1(x1,x2):
    return -6*x1*x2

def df1_dx2(x1,x2):
    return 3*x2**2-3*x1**2

def df2_dx1(x1,x2):
    return 3*x1**2-3*x2**2

def df2_dx2(x1,x2):
    return -6*x1*x2

cnt=0
not_converge=0
x=np.array([-0.5,0.5]) 
origin_x=np.array(x)

while True:
    df = np.array( [ [ df1_dx1( x[0], x[1] ) , df1_dx2( x[0] , x[1]) ] , [ df2_dx1( x[0] , x[1] ), df2_dx2( x[0], x[1] ) ] ] )
    if(df1_dx1(x[0],x[1])==df1_dx2(x[0],x[1])) and (df2_dx1(x[0],x[1])==df2_dx2(x[0],x[1])): #行列のランクが0か１になってしまったと
        print("収束しませんでした")
        plt.scatter(origin_x[0],origin_x[1],c="black")
        break

    f=np.array([f1(x[0],x[1]),f2(x[0],x[1])])
    y=np.linalg.solve(df,-f)
    past_x=np.array(x)
    plt.annotate(not_converge,xy=(x[0],x[1]),size=10)

    for i in range(2):
        x[i]=y[i]+x[i]

    plt.scatter(*x,color="black")
    not_converge+=1

    if not_converge==100: #収束しない
        not_converege=0
        plt.scatter(origin_x[0],origin_x[1],c="black")
        print("収束しませんでした:",*x)
        break

    if abs(x[0]-past_x[0])<=10**-6 and abs(x[1]-past_x[1])<=10**-6:
        if x[0]<0 and x[1]>0:
            plt.scatter(origin_x[0],origin_x[1],c="red")
            print("収束しました(第二象限):初期点(",origin_x[0],",",origin_x[1],")",cnt,"個め:",*x)
        elif x[0]<0 and x[1]<0:
            plt.scatter(origin_x[0],origin_x[1],c="blue")
            print("収束しました(第三象限):初期点(",origin_x[0],",",origin_x[1],")",cnt,"個め:",*x)
        else:
            plt.scatter(origin_x[0],origin_x[1],c="green")
            print("収束しました( x軸上  ):初期点(",origin_x[0],",",origin_x[1],")",cnt,"個め:",*x)

        break

xrange=np.arange(-4,4,0.1)
yrange=np.arange(-4,4,0.1)
x,y=np.meshgrid(xrange,yrange)
plt.axis([-4,4,-4,4])
plt.gca().set_aspect('equal', adjustable='box')

z1=y**3-3*x**2*y
z2=x**3-3*x*y**2-4

plt.contour(x,y,z1,[0],colors="red")
plt.contour(x,y,z2,[0],colors="blue")

plt.show()