ニュートン法の例題で方程式の解と近似値を計算するPython

ニュートン法の例題で方程式の解と近似値を計算するPython

記事内に広告を含む場合があります。
tmp
ニュートン法の例題と解析の要点
📐
接線を利用した近似

曲線の接線を用いて、解（x軸との交点）に素早く近づくアルゴリズムの基本原理。

🏗️
構造解析での応用

非線形解析における荷重変形曲線の追跡に使われ、剛性マトリクスの更新に対応します。

⚠️
収束しないリスク

初期値の選び方によっては解に収束せず、振動や発散を起こす「カオス」な挙動を示します。

ニュートン法の例題

ニュートン法の基本原理と接線の式

 

ニュートン法（ニュートン・ラフソン法）は、方程式 $f(x) = 0$ の解を数値的に求めるための反復計算アルゴリズムです。建築分野における構造解析、特に大変形解析や材料非線形解析などの「非線形問題」を解くエンジンの役割を果たしており、実務で解析ソフトを使う際にも、その挙動を理解しておくことは非常に重要です。
この手法の核心は、「曲線上の点における接線を引き、その接線とx軸との交点を新しい近似解とする」 という幾何学的なアプローチにあります 。
参考）ニュートン法の解説とそれを背景とする入試問題
​
具体的な手順は以下の通りです。


  1. 初期値の設定: 解に近いと思われる適当な値 $x_0$ を定めます。

  2. 接線の計算: 曲線 $y = f(x)$ 上の点 $(x_0, f(x_0))$ における接線を求めます。この接線の傾きは導関数 $f'(x_0)$ で表されます。

  3. 交点の算出: 接線がx軸と交わる点（$y=0$ となる点）を求め、これを新しい近似値 $x_1$ とします。

  4. 反復: この操作を繰り返し、$\Delta x$（修正量）が許容誤差以下になるまで続けます。

これを数式で表すと、以下の漸化式が導き出されます。
xn+1=xn−f(xn)f′(xn)x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}xn+1=xn−f′(xn)f(xn)
この式は、構造力学における「荷重＝剛性×変位 ($F = K \cdot \delta$)」の関係において、剛性 $K$ が接線勾配 $f'(x)$ に、不釣り合い力が $f(x)$ に相当すると考えると、直感的に理解しやすくなります。ニュートン法は「2次収束」という特性を持ち、解に近い初期値を選べば、非常に高速に真の値へ収束するという強力なメリットがあります 。​
非線形の問題は一体どのように解くの？ 解析ソフトはどんな計算をしているのか
参考）https://monoist.itmedia.co.jp/mn/articles/2212/14/news007.html
​

ニュートン法でルート2の近似値を計算する例題

ここでは、電卓もコンピュータも使わずに手計算でニュートン法の威力を体感できる古典的な例題、「$\sqrt{2}$ の近似値を求める」計算を行います。これは $x^2 = 2$、つまり $x^2 - 2 = 0$ という方程式の正の解を求めることと同義です。
設定する関数と導関数は以下のようになります。


  • 対象関数: $f(x) = x^2 - 2$

  • 導関数: $f'(x) = 2x$

漸化式にこれらを代入すると、以下のようになります。
xn+1=xn−xn2−22xn=2xn2−(xn2−2)2xn=xn2+22xn=12(xn+2xn)x_{n+1} = x_n - \frac{x_n^2 - 2}{2x_n} = \frac{2x_n^2 - (x_n^2 - 2)}{2x_n} = \frac{x_n^2 + 2}{2x_n} = \frac{1}{2}\left(x_n + \frac{2}{x_n}\right)xn+1=xn−2xnxn2−2=2xn2xn2−(xn2−2)=2xnxn2+2=21(xn+xn2)
この式を使って、初期値 $x_0 = 2.0$ からスタートして計算してみましょう。



  • 1回目の反復 ($n=0$):


    • $x_0 = 2.0$


    • $x_1 = \frac{1}{2}(2.0 + \frac{2}{2.0}) = \frac{1}{2}(2.0 + 1.0) = 1.5$


    • この時点での2乗値: $1.5^2 = 2.25$（まだ誤差が大きい）




  • 2回目の反復 ($n=1$):


    • $x_1 = 1.5$


    • $x_2 = \frac{1}{2}(1.5 + \frac{2}{1.5}) = \frac{1}{2}(1.5 + 1.3333...) \approx 1.41666...$


    • この時点での2乗値: $1.41666^2 \approx 2.0069$（急速に2に近づいています）




  • 3回目の反復 ($n=2$):


    • $x_2 \approx 1.41666667$


    • $x_3 = \frac{1}{2}(1.41666667 + \frac{2}{1.41666667}) \approx 1.41421568...$


    • 実際の $\sqrt{2} = 1.41421356...$ と比較すると、小数第5位まで一致しています。




たった3回の計算で、実務上ほぼ問題ない精度まで収束しました。これがニュートン法の強みです。もし「二分法」のような単純な探索を行った場合、同じ精度を出すには数十回の反復が必要になることもあります 。
参考）https://www2.akita-nct.ac.jp/saka/Lecturenote/lecture/5e/text/11.pdf

​

ニュートン法の方程式をPythonで実装するコード

次に、このロジックをPythonで実装してみましょう。建設DXや業務効率化の文脈で、独自の数値解析ツールを作成する場合、ライブラリに頼るだけでなく、基礎ロジックを自分で記述できることはデバッグ能力向上につながります。
以下は、任意の方程式に対してニュートン法を適用するPythonコードの例題です。

pythonimport math


 


def newton_method(f, df, x0, tol=1e-8, max_iter=50):
"""
ニュートン法による解の探索


 


Parameters:
f       : 対象とする関数 f(x)
df      : 導関数 f'(x)
x0      : 初期値
tol     : 許容誤差 (tolerance)
max_iter: 最大反復回数


 


Returns:
float: 収束した解
"""
x = x0
for i in range(max_iter):
y = f(x)
y_prime = df(x)


 


# 導関数が0に近い場合、ゼロ除算エラーを防ぐ
if abs(y_prime) < 1e-10:
print("エラー: 接線の傾きが0になりました。収束しません。")
return None


 


# ニュートン法の更新式: x_new = x - f(x) / f'(x)
step = y / y_prime
x_new = x - step


 


print(f"{i+1}回目: x = {x_new:.10f}, 誤差 = {abs(step):.10e}")


 


# 収束判定
if abs(step) < tol:
print("収束しました。")
return x_new


 


x = x_new


 


print("指定された回数内で収束しませんでした。")
return x


 


# 実装のテスト: x^2 - 2 = 0 を解く
if __name__ == "__main__":
# 関数定義: f(x) = x^2 - 2
func = lambda x: x**2 - 2
# 導関数定義: f'(x) = 2x
deriv = lambda x: 2*x


 


# 初期値 2.0 で実行
result = newton_method(func, deriv, 2.0)
print(f"最終結果: {result}")
print(f"検算(result^2): {result**2}")

このコードのポイントは、abs(step) < tol による収束判定と、y_prime が0になった場合の例外処理です 。
参考）Pythonでのニュートン法の実装 - PyDocument
​
実務でScipyなどのライブラリを使用する場合は、scipy.optimize.newton を利用するのが一般的ですが、内部でこのような反復計算が行われていることを理解しておくと、パラメータ設定（tolやmax_iter）の勘所が掴みやすくなります 。
参考）Pythonでやってみた(Engineering)：方程式の…
​
Pythonでのニュートン法の実装手順と解説​

ニュートン法の建築構造計算における収束判定と非線形解析

建築構造計算、特に時刻歴応答解析や保有水平耐力計算などの非線形解析において、ニュートン法（ニュートン・ラフソン法）はまさに心臓部と言えます。ここでは、数学的な「方程式の解」が、構造工学的な意味で何を指しているのかを整理します。
構造解析における方程式は、基本的に「力の釣り合い式」です。
R(u)=Fext−Fint(u)=0R(u) = F_{ext} - F_{int}(u) = 0R(u)=Fext−Fint(u)=0
ここで、$F_{ext}$ は外力、$F_{int}(u)$ は変位 $u$ に依存する内力、$R(u)$ はその差である「不釣り合い力（残差力）」です。
ニュートン法を用いて $R(u)=0$ となる変位 $u$ を探す作業こそが、非線形解析の正体です 。
参考）弧長法とニュートン・ラフソン法とは？基本原理や特徴、適用上の…
​
構造計算におけるニュートン法の重要ポイント:


  • 接線剛性マトリクスの更新:
    数学上の $f'(x)$ は、構造解析では「接線剛性マトリクス $K_T$」に相当します。部材が降伏したり、座屈したりして剛性が低下すると、この勾配が変化します。解析ソフトが計算ステップごとに「剛性マトリクスを更新中」と表示するのは、新たな $x_n$ での傾きを計算し直している状態です 。
    参考）ニュートンラプソン法｜CAE・Ansysの活用推進、解析に関…
    ​

  • 収束判定の厳しさ:
    解析設定における「収束許容値（Tolerance）」は、不釣り合い力がどの程度小さくなれば「釣り合った」とみなすかの閾値です。これを厳しくしすぎると計算が終わらず、緩くしすぎると解析精度が落ち、建物が実際よりも強く（あるいは弱く）評価される危険があります。

  • スナップスルー問題:
    アーチ構造やドームの座屈解析などでは、荷重が増えるにつれて剛性がゼロになり、その後負の剛性（スナップスルー）を示すことがあります。この点（極値点）では接線の傾き $f'(x)$ が0になるため、通常のニュートン法では計算が発散し、エラーで止まってしまいます。
    このようなケースでは、ニュートン法の弱点を補うために**「弧長法（Arc-Length Method）」**という、荷重レベル自体も未知数として扱う高度な手法が併用されます 。​

現場で「計算が流れない（収束しない）」というトラブルに直面した際、それが単なるモデル化ミスなのか、構造的な不安定（接線剛性ゼロ）によるニュートン法の限界なのかを見極める視座を持つことは、構造技術者として一歩抜きん出るために不可欠です。
ニュートンラプソン法｜CAE・Ansysの活用推進​

ニュートン法の初期値設定でカオスや振動が発生する罠

最後に、一般的な教科書ではあまり触れられない、ニュートン法の「闇」の部分、すなわち初期値設定によるカオス的挙動について解説します。
ニュートン法は「局所的」には最強の収束速度を誇りますが、「大域的」には収束が保証されていません。初期値 $x_0$ の選び方ひとつで、解にたどり着くどころか、永久にループしたり、全く無関係な解へ飛んでいったりすることがあります 。
参考）http://www.ep.sci.hokudai.ac.jp/~gfdlab/comptech/y2016/resume/070_newton/2014_0626-mkuriki.pdf
​
特に興味深いのが**「ニュートン・フラクタル（Newton Fractal）」**と呼ばれる現象です。
例えば、複素平面上で $z^3 - 1 = 0$ のような単純な方程式を解く際、どの初期値が3つの解（$1, e^{i2\pi/3}, e^{i4\pi/3}$）のどれに収束するかを色分けしてプロットすると、境界部分に驚くほど複雑な自己相似形のフラクタル模様が現れます。
これは実数の方程式でも起こり得ます。以下のようなケースで発生します。


  • 変曲点近傍の罠:
    変曲点（グラフの凹凸が変わる点）付近を初期値に選ぶと、接線が遠くへ飛びすぎてしまい、解から遠ざかる発散挙動を起こします。

  • 循環（リミットサイクル）:
    $x_0$ での接線が $x_1$ を指し、$x_1$ での接線がまた $x_0$ に戻るような特殊な関数形状の場合、永遠にその2点間を振動し続け、無限ループに陥ります。構造解析で言えば、荷重ステップの途中で変位が行ったり来たりして収束しない現象がこれに近い状態です 。
    参考）http://www.yamamo10.jp/yamamoto/lecture/2007/5E_comp_app/nonlinear_equation/nonlinear_eq_html/node4.html
    ​

  • 導関数が0になる点:
    前述の通り、$f'(x)=0$ となる極値付近では、接線が水平になり、次の解の候補が無限遠点となって計算不能（Overflow）に陥ります。

対策としての減衰ニュートン法:
このような振動や発散を防ぐため、実務的なコードでは更新量を $\Delta x = -f(x)/f'(x)$ そのまま使うのではなく、$\Delta x = -\alpha (f(x)/f'(x))$ のように、ステップ幅を調整する係数 $\alpha (0 < \alpha \le 1)$ を導入する「減衰ニュートン法」が使われることがあります。
「ライブラリを使えば答えが出る」と過信せず、解こうとしている方程式の形状（構造特性）によっては、初期値の変更やステップ幅の調整といった泥臭い試行錯誤が必要になることを、プロフェッショナルとして心に留めておく必要があります。
非線形方程式の数値計算と収束しない例​

 

 


電力: 現状と新発電法 (ニュートンムック Newton別冊)