【プログラミング言語Julia】パッケージ［Symbolics.jl］で数式処理を組み合わせたニュートン法を実装する。

今回もJuliaの数式処理パッケージSymbolics.jlの事例紹介になります。過去の記事は以下になります。

otepipi.hatenablog.com

また、記事内容は以下書籍のプログラムを参考にしています。

MATLABとOctaveによる科学技術計算

丸善出版

Amazon

ニュートン法とは

ニュートン法は $f(x)=0$ の解 $x$ を求めるための著名な方法の一つで、以下の繰り返し計算によって求める事ができます。

\begin{equation} x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} \end{equation}

今回はJuliaのパッケージSymbolics.jlを使ってニュートン法を実装してみます。

Symbolics.jlを使うことの利点

$f(x)$ を手動でプログラムに入力すれば、数式処理システムにより数式 $f'(x)$ が自動で求まるからです。

数式処理システムを使用しない場合、人手で $f(x)$ を微分して $f'(x)$ の数式を求めて入力することになります。人手による計算を可能な限り少なくすることは重要です。また近年では数値計算と数式処理の組み合わせによるシミュレーションが増えており、今回はその流れに乗って数値計算のニュートン法に数式処理を活用する試みになります。

パッケージの宣言

Symbolics.jlを宣言します

using Symbolics

数式の定義

今回、練習として以下の数式の解を求めることにします。

\begin{equation} f(x) = x + e^{x} + \frac{10}{1 + x^{2}} - 5 \end{equation}

@variables x
f = x + exp(x) + 10/(1+x^2) - 5

\begin{equation} -5 + x + e^{x} + \frac{10}{\left( 1 + x^{2} \right)} \end{equation}

数式の微分を求める。

Differential(x)を利用して数式の微分を求めます。expand_derivativesを用いて微分した数式を展開しています。

D = Differential(x);
df = expand_derivatives(D(f))

\begin{equation} 1 - \frac{20 x}{\left( 1 + x^{2} \right)^{2}} + e^{x} \end{equation}

$f,df$ を使ってニュートン法を使用する。

$diff = \frac{f(x^{(k)})}{f'(x^{(k)})}$ として先に求めておきます

diff = f/df

\begin{equation} \frac{\left( -5 + x + e^{x} + \frac{10}{\left( 1 + x^{2} \right)} \right)}{\left( 1 - \frac{20 x}{\left( 1 + x^{2} \right)^{2}} + e^{x} \right)} \end{equation}

今回初期値x0=0.3、最大の繰り返し数maxniter = 1000, 許容誤差tol = 0.0001としてみます。最大繰り返し数に達したにも関わらず許容誤差が大きい場合はエラーとします

x0 = 0.3
maxniter = 1000
tol = 0.0001

niter = 0;

xk = x0
err = tol + 1

while abs(err) > tol && niter < maxniter
    niter = niter + 1; 
    xk =xk - substitute(diff,Dict([x => xk]))
    err = substitute(f,Dict([x => xk]))
end

abs(err) > tol && error("over max iteration") ;