非線形方程式の近似根を求めるライブラリに頼ってみよう

自作の Newton 法ルーチンがうまくうごかないような問題もあるだろう．それに、Newton 法は微分関数を計算して求めておく必要がある．これはときに大変に面倒だし、そしてバグのもとにもなりかねない．

  • 数式処理ソフトなどを使うことも出来るし、「自動微分」という良い方法もあるけど、そうした方法についてはまた別のところで触れよう．

そこで、既存のライブラリを使ってみよう．

今回は、NLsolve という、非線形方程式の近似根を求めるそのものズバリのライブラリ(Package)を使ってみよう．

パッケージの使い方

julia には, 公開・配布されているライブラリがたくさんあり、これらはパッケージ(package)と呼ばれている．

julia package list

にいくと、オフィシャルなリストが見られるので、暇な時に見ておくと良いだろう(ちなみに、現時点でパッケージの数は 1500 を超えている)．

こうしたパッケージの使い方はシンプルで、インストールして、あとは使うたびに呼び出すというだけだ． ここでは、NLsolve パッケージを例にして説明しよう．

インストール (一回だけやれば良い)

julia 上で Pkg.add 命令を使ってインストールできる．この場合はパッケージ名 NLsolve を与えて

Pkg.add("NLsolve")

とすることになる(今、やろう!)．

  • JuliaBox "Next" からはパッケージは GUI で管理できるようになっている．管理画面の上にある "Packages" をクリックしてみよう！

パッケージの呼び出し

パッケージの呼び出しは他の言語のライブラリと同様で、そのパッケージの機能を使う前にプログラム中で一度「使うよ」と宣言しておけば良い． julia の場合は、

using パッケージ名

として呼び出すことになる．というわけで、NLsolve のインストールが済んでいる人は、次のようにして使用宣言を今しておこう．

In [1]:
using NLsolve

NLsolve パッケージの使い方

このパッケージは非線形問題の汎用解法ライブラリなので、そのままだとちょっと使いにくい． 詳しくは

NLsolve: 非線形方程式の求根 package

などを見てもらうとして、ここでは処方箋だけ書いておこう．

この処方箋とは、using NLsolve のあとに、続けて nls という関数を julia で定義しておいて、あとはこの関数 nls を使うという手である． このためには、この関数を定義しておかないといけないので、やっておこう．次のようになる．

In [2]:
function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vin,vout)->vout[1]=func(vin[1],params...), [ini])
        v = r.zero[1]
    else
        r = nlsolve((vin,vout)->vout .= func(vin,params...), ini)
        v = r.zero
    end
    return v, r.f_converged
end
Out[2]:
nls (generic function with 1 method)

関数 nls の使い方

スカラー問題と、ベクトル問題で分けて解説しておこう．

スカラー問題の場合

解きたい方程式が, スカラー変数 $x$ に対してスカラー値関数 $f$ を用いて $f(x,params)=0$ ($params$ はパラメータ, つまり、関数 $f$ の２つ目以降の引数) となっている場合は、

nls( f, params, ini = xの初期値 )

とすれば近似根を取り出せる． 文法的に、初期値設定 ini = は省略できないが、params はあってもなくても複数でも良い．

ベクトル問題の場合

解きたい方程式が, ベクトル値変数 $\boldsymbol{x}$ に対してベクトル値関数 $\boldsymbol{f}$ を用いて $\boldsymbol{f}(\boldsymbol{x}, params) = \boldsymbol{0}$ (params はパラメータ、つまり、関数 $f$ の２つ目以降の引数) となっている場合は、

nls( f, params, ini = xの初期ベクトル )

とすれば近似根ベクトルを取り出せる． 文法的に、初期値設定 ini = は省略できないが、params はあってもなくても複数でも良い．

スカラー問題をやってみる

パラメータがない、とても簡単な例をやってみよう． $\cos$ 関数の根、つまり、$\cos(x) = 0$ を満たす $x$ を一つ求めてみよう． それには、例えば近似根の初期値を $x_0 = 1.0$ などとして次のようにすれば良い．

In [3]:
nls(cos, ini = 1.0)  # ini = 1.0 は近似根 x の探索初期値．
Out[3]:
(1.570796326795489, true)

すると上のように出力され、一つ目の数字が近似根で、２つめの "true" が「この近似根で $f \cong 0$ となっている」ことが確かめられる．

  • 初めて nls を動かす時は結構待たされる．JIT という機構のためだ．まあそういうもんだと思っておこう．

一応、近似根がうまく求まっていることを確認しておこう．

In [5]:
x = nls(cos,ini = 1.0)[1]  # 並んでいるデータの一つ目だけ取り出す．"true" は不要なので．
cos(x)  # そして、cos に代入して、ゼロに近いかどうかを見てみる．
Out[5]:
-5.923537736000261e-13

結果は 0 に大変近い．確かに、この近似根は十分「良い」と言って良さそうだ．

次に、パラメータ(関数 $f$ の２つ目以降の引数のこと)が有る例をやってみよう． 例として、パラメータ $p$ の平方根を求めるような問題を考えてみよう． それには､ 関数を $f(x,p) = x^2 - p$ としてこの根 $x$ を求めればいい． だから、次のようになるだろう．

In [6]:
f(x,p) = x^2 - p

nls(f, 2.0, ini=1.0) # 二つ目の引数 2.0 はパラメータ．ini = 1.0 は近似根 x の探索初期値．
Out[6]:
(1.4142135623746899, true)

すると確かに $\sqrt{2}$ 相当が得られている．

パラメータは複数でも良いので、それを例でみてみよう． 例えば、上の関数を少しだけいじって $f(x,p,q) = x^q - p$ としてこの根 $x$ を求める問題を考えよう．

In [7]:
f(x,p,q) = x^q - p

nls(f, 2.0, 2.1, ini=1.0) # x^(2.1) - 2.0 = 0 の近似根を求める
Out[7]:
(1.391065619249075, true)

この値が確かに近似根になっていることは、結果の "true" 部分からも分かるが、確かめてみると、

In [9]:
x = nls(f, 2.0, 2.1, ini=1.0)[1]

x^(2.1) - 2.0
Out[9]:
9.798828415341632e-12

となるので、確認できる．

近似根がうまく求まらない問題も扱ってみよう． 単純なケースとして、負の数字の平方根を求めることに相当する近似根計算をさせてみよう．

In [12]:
f(x,p) = x^2 - p

nls(f, -0.0001, ini = 1.0)  # x^2 = - 0.0001 を解こうとさせている．当然、実数解は存在しないので…
Out[12]:
(-4.385150012659097e-11, false)

出力をみると、(探索ルーチンが頑張って求めた)一見近似根っぽい数字が一つ目の出力項に出ているけれども、2つめの出力項の false によって、それは $f \cong 0$ を満足するような近似根ではない、という情報がわかる．つまり、

まともな近似解は求まらなかったよ！ と言っている のである．

「ややこしい関数で近似根を求めさせると失敗するかも」というようなケースでは、この「計算に成功したかどうか」をきちんとチェックするようにプログラムを書いたほうが良さそうだ．

ベクトル値問題もやってみよう

入出力値がベクトルの問題もやってみよう． 例として，パラメータ $p,q$ に対して，

$\left\{\begin{array}{rcl} \cos^2(x) + \sin^2(y) & = & p, \cr \sin^2(x) - \cos^2(y) & = & q \end{array}\right.$

という連立方程式を満たす近似根 $x,y$ が欲しいとしよう．

まず，この問題に対して等価な $h(\boldsymbol{x}, p, q) = \boldsymbol{0}$ という式を満たすように関数 $h$ を次のように定義しよう．

In [13]:
function h(u, p, q) 
    
    a,b = u[1],u[2] # u[1] とか毎回書くのは面倒なので，a,b で代用できるようにして…
    
    return [ cos(a)^2 + sin(b)^2 - p, sin(a)^2 - cos(b)^2 - q ]
    
end
Out[13]:
h (generic function with 1 method)

試しに，適当な値を入れてチェックしてみる．

In [14]:
h([4.0,4.0], 0.5, 0.4)
Out[14]:
2-element Array{Float64,1}:
  0.5   
 -0.2545

確かに、ベクトル値が出力され、動作はきちんとしていそうだ．

では，上のご託宣に従って，近似根を求めてみよう．

In [15]:
p = 0.5
q = 0.4
ini_v = [4.0,4.0]

z = nls(h, p, q, ini = ini_v) # パラメータが複数でも大丈夫．
Out[15]:
([4.48688, 3.87691], true)

true と出ているので大丈夫なはずだが、この値が近似根になっていることを確認しておこう．

In [16]:
h(z[1], p, q)  # 結果を格納した z の一つ目が、近似解ベクトルなので [1] をつけるよ．
Out[16]:
2-element Array{Float64,1}:
  1.36335e-13
 -1.36668e-13

確かにほぼ ゼロベクトルが出力されるので、OK と言って良さそうだ．

自力でやってみよう

$sin(x) = 3/4$ の近似解を求めてみよう

ほぼ即座にできるだろう．そこで、できれば、

  • 上の例を見ないでやってみよう
  • Newton 法の結果と比較してみよう

というチャレンジをしてみよう．

非線形連立方程式の近似解も求めてみよう

これも Newton 法のときにでてきた問題だが、

$\left\{\begin{array}{rcl} x + 2y & = & 0, \cr x^2 + 2y^2 & = & 0 \end{array}\right.$

の近似解を求めてみよう．そして、やはり

  • 上の例を見ないでやってみよう
  • Newton 法の結果と比較してみよう

というチャレンジをしてみよう．