任意の方程式（関数）に対してニュートン法で1つの解を求めるプログラムを1から書いてみた。
欠陥あるかもしれません。

ニュートン法

$x_{new}=x_{old}-\frac{f(x_{old})}{f(x_{old})^{\prime}}$

設定するもの

・解きたい方程式（関数）
・初期値（真の解に近めが望ましいが、そんなんに気にしない。）
・許容誤差

def newton_method(func, initial_x, r_error):

こういう書き出しで作ってみます。

更新の流れ

$x_{new}=x_{old}-\frac{f(x_{old})}{f(x_{old})^{\prime}}$

に基づいて、

$|f(x_{初期値})| >$ 許容誤差 　ならば

$x_{初期値}-\frac{f(x_{初期値})}{f(x_{初期値})^{\prime}}=x_{更新値1}$

$|f(x_{更新値1})| >$ 許容誤差 　ならば

$x_{初期値}$ の部分を $x_{更新値1}$と改めた上で　次の 更新値2 を求める。

$x_{更新値1}-\frac{f(x_{更新値1})}{f(x_{更新値1})^{\prime}}=x_{更新値2}$

繰り返していけばそのうち

$|f(x_{更新値n})| <$ 許容誤差 となり許容誤差の範囲に入る。

関数次第では収束しない場合もある。

実装1

# 実装
def newton_method(func, initial_x, r_error):
    """
    funcの部分はlambda関数で設定する。
    func = lambda x : x**2 - 7
    """
    count = 0
    while True:
        count += 1
        
        y_f = func(initial_x)

        # 許容誤差との比較
        if abs(y_f) < r_error:
            break
            
        # 許容誤差を満たさない場合は更新
        # 微分係数を求めて更新
        
        # 微分係数は中心差分を用いる
        y_d = (func(initial_x + 1e-10) - func(initial_x - 1e-10)) / (2 * 1e-10)
        
        # 更新値を求める
        x_new = initial_x - (y_f / y_d)
        print("{}回目:X_new = {}".format(count,x_new))
        # 更新
        initial_x = x_new

関数次第で無限ループします。
実装２で離脱する設定を追加しています。

試運転

$\sqrt{7}の近似$
方程式（関数）:$x^2-7=0$
初期値は適当に30
許容誤差:0.001

# 関数部分はlambda関数で指定する
newton_method(lambda x : x**2 - 7, 30, 0.001)

実装の補足

代入処理

$f(x_{old})$の代入処理について。

今回は関数が最初から決まっておらず、任意に決定できるようにしてあります。

このような場合に代入処理は、あんまり直球すぎると処理されません。

余興で、設定した関数に対して、xを代入する、『代入machine』を作ってみます。

def dainyuu_machine(func, x):
    y = func
    return y

dainyuu_machine(x**2-4,3)

これは x is not defined のエラーが返ってきます。

funcの部分は式を突っ込むのではなく、関数を使うことにします。（lambda関数を使います。）

def dainyuu_machine(func, x):
    y = func(x)
    return y

dainyuu_machine(lambda x : x**2-4,3)

これだと $f(3)=3^2-4$ が計算されます。
任意の式を設定してから代入への流れはこの仕組みを使ってます。

lambda関数

def func(a):
    return a + 1
func(3)


# lambda関数
x = lambda a: a + 1
x(3)

微分係数の近似について

鉛筆を持って、$f(x)=x^2-7$を微分するのは造作もないことですが、
プログラムで動かすとなると、何次式になるかなどの場合分けがあり、
サクッと式の形に落として、プログラムとして実装することがなかなか難しそうです。

今回は代入した値が求まればそれで良いので、
導関数の定義式？的なので近似します。
$$f(a)^{\prime}= \frac{f(a+h)-f(a)}{a+h-a}= \frac{f(a+h)-f(a)}{h}$$

これを使っても良いのですが、より誤差が少ない中心差分公式を用いて微分係数を返してみました。
$$f(a)^{\prime}= \frac{f(a+h)-(a-h)}{2h}$$

参考：https://gihyo.jp/dev/serial/01/java-calculation/0070

$a$の部分は$x_{old}$に該当。代入の処理はlambda関数を使用。
hの部分は微小誤差として、勝手に1e-10としていますが、
計算効率に関わってくる部分っぽいので、本当はあんまり勝手に決めない方が良いのかもしれません。

実装2

・countと誤差の推移のグラフを追加
・100回で収束しない場合の撤退を追加

import matplotlib.pyplot as plt
import numpy as np


# 実装
def newton_method(func, initial_x, r_error):
    """
    funcの部分はlambda関数で設定する。
    func = lambda x : x**2 - 7
    """
    count = 0
    # 折線表示のために誤差結果の保存
    y_f_value_list = []
    
    while True:
        count += 1

        y_f = func(initial_x)
        # 誤差リストに追加
        y_f_value_list.append(abs(y_f))
        
        # 許容誤差との比較
        if abs(y_f) < r_error:
            break
        # 許容誤差を満たさない場合は更新
        # 微分係数を求めて、更新
        
        
        # 微分係数は中心差分を用いる
        y_d = (func(initial_x + 1e-10) - func(initial_x - 1e-10)) / (2 * 1e-10)

        # 更新値を求める
        x_new = initial_x - (y_f / y_d)
        print("{}回目:X_new = {}".format(count,x_new))
        # 更新
        initial_x = x_new
        
        # 100回tiralで収束しない場合、計算を止める
        if count == 100:
            print("近似解見つからず。")
            break
    
    # countと誤差の推移
    x = np.arange(1, count+1, 1) # 1 ~ n回
    y = np.array(y_f_value_list) # 1 ~ n回の各誤差
    plt.plot(x, y)
    # 横軸と縦軸のラベルを追加
    plt.xlabel('trial')
    plt.ylabel('r_error')
    
    # 許容誤差の横線
    plt.hlines(r_error, 1, count,color="red",linestyles='dashed')
    # 折れ線と許容誤差のボーダーを一括で表示
    plt.show()

$x^2-7=0$ を初期値20から近似する。

newton_method(lambda x : x**2 - 7, 20, 0.001)

実装3

実装1,2は許容誤差を高さに設定して、近似している。

ここではニュートン法による横の動きの更新が浅くなったら近似解にする感じのプログラムも書いてみた。

今回、関数は最初に外に出して定義。
x_ini:初期値
JJ:繰り返し回数
EPS:許容誤差

$\frac{|x_{new}-x_{old}|}{|x_{old}|} > EPS $

で評価する。

$|x_{new}-x_{old}|$が小さくなればなるほど収束していると考えられる。
値が極端に大きい時にも対応できるように$|x_{old}|$で正規化したものを使った。

def ff(x):
    return x**2-7

# 微分係数を返す関数（中心差分を用いる）
def df(x):
    return (ff(x + 1e-10) - ff(x - 1e-10)) / (2 * 1e-10)

def newton_method(x_ini, JJ, EPS):
    x_old = x_ini
    for i in range(JJ):
        x_new = x_old - ff(x_old)/df(x_old)
        # 評価部分
        if (abs(x_new-x_old)/abs(x_old)) < EPS:
            break
        
        else:
            x_old = x_new
            
    if i == JJ-1:
        print("Do not CONVERGE")
    else:
        return x_new
    
if (__name__=='__main__'):
    JJ = 100
    EPS = 0.001
    x_ini = 20
    x_sol = newton_method(x_ini, JJ , EPS)
    print("solution:{}".format(x_sol))