gnuplotで方程式を解く

方程式を解くための初歩的な手法であるNewton法は、gnuplotでも使うことができます。

Newton法の解説

Newton法とは、方程式を数値的に解くためのもっとも簡単な手法の一つで、以下のように「幾何学的」に解を求めるものです。

ここでは、ある関数f(x)についての以下のような方程式を考えます。

まず、xの初期値x0を（適当に）定めます。整数n = 0, 1, 2, ...に対し、x=xnでのfの値f(x0)と、x=xnでのfの微分f'(xn)=df/dx(xn)を求め、xn+1を以下の式で定めます。

これは、下図のように、x=xnにおける接線がy=0の直線を横切る点をxn+1と定めていることになります。新しいxn+1に対して同様の操作を行い、xn+2を求め、さらに繰り返していくと、初期値と関数が適切であれば、数列xnは上記の方程式の解xsolに収束していきます。

Newton法の実装その1（ループを使う方法）

まずは、gnuplotのループ機能を使ってNewton法による解探索を実行する方法を解説します。以下では、log(x)-1という関数のゼロ点を求めることを考えます。方程式log(x)-1=0の解は、x = e = 2.71828...です。whileループの中の青字で示した部分で、先に説明した次の点を求める方法を用いています。

f(x)=log(x)-1 # 解（ゼロ点）を求める関数

dx =  1E-6       # 微分をする時に使う微小数
EPS = 1E-8       # 解の収束判定条件
LOOP_LIM = 1000  # ループ回数の上限
x0=1.0           # 探索の初期値

set xrange [0.5:5]
set yrange [-2:2]
set xzeroaxis

# Newton法ループ
i=0
while (abs(f(x0)) > EPS && i < LOOP_LIM){
     df = (f(x0+dx)-f(x0-dx))/(2*dx)
     set title sprintf("loop: %d, x0 = %g, f(x0) = %g", i, x0, f(x0))
     set arrow 1 from first x0,f(x0)+0.5 to first x0,f(x0)+0.05 size screen 0.02,15 filled lw 2
     plot f(x), df*(x-x0)+f(x0)
     pause -1
     if (abs(f(x0)) > EPS) {x0 = -f(x0)/df + x0}
     i=i+1
     }

# 最終的に得られた解の表示
set title sprintf("Final result (after %d loops): \nnumerical solution: x = %.10f\n", i-1, x0, f(x0))
set arrow 1 from first x0,f(x0)+0.5 to first x0,f(x0)+0.05
plot f(x), df*(x-x0)+f(x0)

print "numerical solution:  x = ", x0
print "exact solution: exp(2) = ", exp(1)
print "error: ", x0 - exp(1)

この例を実行すると、以下のようにNewton法の各ステップの様子をグラフで見ながら、xnがに向かって収束していく様子を見ることができます。

最後に以下のように得られた解を表示して終了します。

実際に得られた解が、解析的な解x = eとよく一致していることがお分かりいただけると思います。

Newton法の実装その2（再帰定義を使う方法）

再帰定義をうまく工夫することでも、Newton法を利用することができます。以下も、先の例と同じくlog(x)-1=0の解を求めていますが、再帰定義を使っています。newton_f関数が、f(x)=0の解を求めるための関数です。この定義の部分では、関数f(x)の絶対値がEPSよりも大きければ-f(x)/((f(x+dx)-f(x-dx))/(2*dx)) + xを新しい変数にして再度newton_f関数を呼び出しています。もし小さければ、x自身を解として返します。

f(x) = log(x)-1 # 解（ゼロ点）を求める関数

dx =  1E-6       # 微分をする時に使う微小数
EPS = 1E-8       # 解の収束判定条件
x0=1.0           # 探索の初期値

# 再帰定義による解
newton_f(x) = (abs(f(x)) > EPS ? newton_f(-f(x)/((f(x+dx)-f(x-dx))/(2*dx)) + x) : x)

# 解を表示する
print "solution (newton): x = ", newton_f(x0)
print "solution (theory): x = ", exp(1)

# wgnuplotを表示させるために必要。ほかのターミナルならば不要。
plot x

このスクリプトを実行すると、以下のように解を表示します。確かに、変数EPSで指定した精度で解が求まっています。

この方法の面白いところは、解を求める部分が関数になっている点です。したがって、これをグラフのプロットなどにも使うことができますが、そのためにはfやnewton_fの引数の数を増やして、プロットに使うための変数を入力できるようにしなければなりません。

例えば、以下の例では、log(x)-a=0の解をaの関数としてプロットしています。これは当然exp(a)と一致するはずです。実行例を見てみると、確かにnewton_f(x0,x)はexp(x)と一致していることがわかります。

f(x, a) = log(x)-a

EPS = 1E-8
dx = 1E-6
x0 = 1.0

# 再帰定義による解
newton_f(x,a) = (abs(f(x,a)) > EPS ? newton_f(-f(x,a)/((f(x+dx,a)-f(x-dx,a))/(2*dx)) + x, a) : x)

set xrange [1:5]
plot exp(x), newton_f(x0,x)