[Python]多変数非線形方程式をNewton-Raphson法で解く

Python

\begin{align} f(y,x)&=\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}

をNewton-Raphson法を使って解いてみましょう。

処理の流れ

$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を更新する$

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

実装してみよう

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


x=np.array([random.uniform(-1,1),random.uniform(-1,1)])
print("初期点: ",*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("収束しませんでした")
    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:
    print("収束しました",*x)
    break

"""実行結果
初期点: -0.8994512575217808 -0.5665343822593587
収束しました -0.7937005259840998 -1.3747296369986026

初期点: 0.599898843152864 0.22199829203012378
収束しました 1.5874010519685833 -4.74652956283883e-13

初期点: -0.9391760864852678 0.8686860309766029
収束しました -0.7937005259842101 1.3747296369986743
"""

収束点は3つあると思います。

もちろんですが、収束する点は方程式の交点なので、数などは方程式によって違います。

この内容はある大学の課題だったそうですね笑

もしこの記事をみてくれいる人で役立ったと思ったら是非コメント、いいねをしてくれると嬉しいですね!

関連記事も載せておきますね!

今回は以上になります。

ではまた。

コメント

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