



実行結果はこんな感じになりました。
実装してみよう
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か1になってしまったと
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()
コメント