[Python] Newton-Raphson法 初期点と収束先の関係性 連立方程式を解く

Python

\begin{align} f(x_{1},x_{2})&=\left\{\begin{array}{ll} x_{2}^3 &- 3x_{1}^2x_{2} \\ x_{1}^3 &- 3x_{1}x_{2}^2 &- 4\end{array}\right.\end{align}

の初期点と収束先の関係性を調べてみましょう!!!

復習

$x^{(k+1)}-x^{(k)}=-f'(x^{(k)})^{-1}f(x^{(k)})、k=0,1,2… $

$変数変換y^{(k+1)}=x^{(k+1)}-x^{(k)}とすると$

$両辺に左からf'(x^{(k)})をかけて$

$連立一次方程式 f'(x^{(k)})y^{(k+1)}=-f(x^{(k)})となる$

上記の方程式を解いて

$y^{(k+1)}=x^{(k+1)}-x^{(k)}$

$x^{(k+1)}=x^{(k)}+y^{(k+1)}としてxを更新する$

こうすることで逆行列の計算をしなくてすみ楽に計算できます。

$f'(x)はヤコビ行列になっています。$

どこが収束先か?

上記の写真のように収束先は3種類あり、方程式の交点となっています。

実装してみる

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
fig=plt.figure()
ax=fig.add_subplot(1, 1, 1)
not_converge=0

while True:

    x=np.array([random.uniform(-1,1),random.uniform(-1,1)])
    #x=np.array([-0.5,0.5])
    origin_x=np.array(x)
    cnt+=1
    while True:
        not_converge+=1
        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("収束しませんでした")
            not_converge=0
            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)
        for i in range(2):
            x[i]=y[i]+x[i]
        
        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")   
            elif x[0]<0 and x[1]<0:
                plt.scatter(origin_x[0],origin_x[1],c="blue")    
            else:
                plt.scatter(origin_x[0],origin_x[1],c="green")
            not_converge=0
            break

    if cnt==10**3:
        break

plt.show()

初期点は$-1<x<1$,$-1<y<1$の範囲に限定しています。

実行結果

これは初期点を100万個プロットしたものになります。

これは1万点をプロットしたものになります。

多ければ綺麗なプロットされていますね。

赤が第二象限、青が第三象限、緑がx軸上に収束することを示しています。

皆さんもいろいろな関数で試してみてはどうでしょうか笑

今回は以上となります。

ではまた。

コメント

タイトルとURLをコピーしました