ニュートン・ラフソン法とは
方程式 f(x)=0 の解を、現在の推定点における接線で線形化して次の推定点を決める反復解法。CAE の非線形ソルバーの中核として日常的に走っている、最も古典的で最も実用的な数値手法の一つだよ。
🙋
このシミュレーター、初期値 x_0 を 2 にしたら3反復で根が見つかったのに、x_0 を −2 にすると急に発散しました。同じ方程式なのに何が起きてるんですか?
🎓
ニュートン法は「接線で根を当てに行く」手法だから、接線の傾き f'(x) が小さい場所に来ると、接線が水平に近くなって x 軸と交わる点が遠くに飛ぶんだ。f(x)=x³−2x−5 の f'(x)=3x²−2 は x=±√(2/3)≈±0.816 でゼロになる。x_0=−2 から出発すると、その不安定域を横切る軌道になって、グラフの右端や左端まで一気にジャンプしちゃう。シミュレーターで初期値スライダーを動かしながら下の誤差曲線(log10)を見ると、「滑らかに減るケース」と「振動して上に跳ねるケース」がはっきり区別できるよ。
🙋
なるほど!じゃあ緩和係数 ω のスライダーを 0.5 にしたら、初期値が悪いときでも収束しますか?
🎓
大正解!更新量を ω 倍に抑えると、暴走しがちな反復が落ち着いて収束しやすくなる。これを「ダンプドニュートン法」と呼んで、Abaqus や ANSYS の非線形構造解析でも内部で同じことをやってる。ただし ω を小さくしすぎると、本来の2次収束(誤差が反復ごとに約2乗で縮む)が壊れて1次収束(線形)に劣化する。シミュレーターの誤差曲線で ω=1.0 と ω=0.3 を比べると、傾きの違い(凸に減るか直線的に減るか)がはっきり見えるはず。
🙋
許容差を 1e^−12 にすると反復回数が増えますか?というか、そんなに小さくする意味あります?
🎓
良い質問。純ニュートン法の2次収束は強烈で、1e^−6 から 1e^−12 まで縮めるのに必要な反復は1〜2回程度しか増えない。だから「許容差を厳しく」しても計算コストはほぼ変わらない。ただ実務では、CAE の倍精度浮動小数の限界(およそ 1e^−15)と、残差の評価誤差(数値積分・モデル不確かさ)の方が先にボトルネックになる。だから線形構造解析なら 1e^−6、非線形なら 1e^−4 〜 1e^−5 を相対残差で取るのが普通だよ。
物理モデルと主要な数式
ニュートン・ラフソン法は、滑らかな関数 f(x) を現在の推定点 x_n まわりで1次テイラー展開し、その線形近似がゼロになる点を次の推定値とする手法。
$$f(x) \approx f(x_n) + f'(x_n)(x - x_n) = 0 \;\Rightarrow\; x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$
緩和係数 ω を導入したダンプドニュートン形式:
$$x_{n+1} = x_n - \omega\,\frac{f(x_n)}{f'(x_n)}, \quad 0 < \omega \le 1$$
2次収束性(誤差 e_n = x_n − x* について):
$$|e_{n+1}| \le \frac{|f''(x^{*})|}{2|f'(x^{*})|}\,|e_{n}|^{2}$$
つまり、根の近くでは反復ごとに有効桁数がほぼ倍増する。これがニュートン法の最大の魅力。
実世界での応用
非線形構造解析(CAE):材料非線形(弾塑性)・幾何非線形(大変形)・接触解析の平衡反復で、残差力 R(u)=F_int(u)−F_ext=0 をニュートン・ラフソン法で解く。剛性行列 K=∂R/∂u が「f'(x)」に相当する。
電気回路シミュレーション(SPICE):非線形素子(ダイオード・MOSFET)を含む直流動作点解析では、Kirchhoff の式 g(V)=0 をニュートン反復で解く。各反復で生成される線形化系(ヤコビ行列)の規模で計算時間が決まる。
機械学習・最適化:勾配 ∇f(x)=0 をニュートン法で解くと、二階導関数 ∇²f(ヘッセ行列)を用いた高速収束が得られる。準ニュートン法(BFGS)はそのヘッセを近似する派生手法。
天体力学・軌道計算:Kepler の方程式 M = E − e·sin E を離心近点角 E について解くとき、ニュートン・ラフソン法が標準的に使われる。衛星追跡や軌道予測の計算で必須。
よくある誤解と注意点
まず一番怖いのは 「2次収束だから初期値はどこでもいい」 という誤解。ニュートン法の2次収束は「根の十分近く」での話で、初期値が悪いと簡単に発散したり別の根にジャンプしたりする。実務では、まず二分法や黄金分割探索で根の存在区間を絞り、そこからニュートン法に切り替える「ハイブリッド法」が安全。Brent 法はそのアイデアを洗練した代表例だ。
次に f'(x) ≒ 0 の場所。接線がほぼ水平になると、更新量 f/f' が異常に大きくなる。シミュレーターで x_0 を 0.8 付近に置くと、f'(0.8)=3(0.64)−2=−0.08 と非常に小さく、次のステップで一気に遠くに飛ぶ様子が見える。これを避けるために、実装ではしばしば |f'| が閾値以下なら ω を縮める、あるいは更新量に上限を設けるトリック(line search、trust region)が組み合わされる。
最後に 「収束した=正解」とは限らない 落とし穴。f(x)=0 が複数の根を持つ場合(このツールの f(x)=x³−2x−5 は実根は1つだけだが、一般の3次式には3つあり得る)、初期値によってどの根に収束するかが変わる。CAE の非線形構造解析でも同じで、座屈・スナップスルー後の解は「正しい平衡解」が複数あり、アーク長法(Riks 法)など特殊な反復で経路を追跡する必要がある。
よくある質問
f'(x) を中心差分 (f(x+h)−f(x−h))/(2h) などの数値微分で置き換える「割線法(Secant method)」や、過去の評価値からヤコビを近似する「ブロイデン法(Broyden)」を使います。CAE では計算コストの大きいヤコビを更新せず使い回す「修正ニュートン法」も一般的です。
重根 x* では f'(x*)=0 となり、純ニュートン法の収束次数は2次から1次(線形)に落ちます。修正式 x_{n+1}=x_n − m·f/f'(m は重複度)を使うと2次収束を回復できます。重複度が事前に分からなければ、Schroder の反復 x_{n+1}=x_n − f·f'/((f')²−f·f'') を使う方法もあります。
はい。ベクトル方程式 F(x)=0 に対しては Δx = −J(x_n)⁻¹ F(x_n) で更新します(J はヤコビ行列)。CAE の非線形有限要素解析がまさにこの形で、J が剛性行列 K に相当します。実装ではヤコビの直接逆行列ではなく、線形ソルバー(LU 分解・反復法)で K·Δu=−R を解きます。
倍精度浮動小数の機械イプシロン(約 2.2e−16)より小さくは原理的に下げられません。さらに、f(x) の評価自体に丸め誤差や数値積分の打切誤差があるため、それを下回る許容差は意味を持ちません。実務では「相対残差 < 1e−6 〜 1e−4」が一般的で、絶対残差ではなく初期残差で正規化して評価します。