- 2018.05 NLsolve package の仕様変更に伴う修正. 以前のコードが動かなくなっていた人は、下記のものに直せば動くはず.
- 2016.12 もう少し簡単に使えるように追記.
非線形方程式の近似根を求めるライブラリに頼ってみよう
Newton 法は、微分関数を「自分が手で計算する必要がある」というデメリットが有る.数式処理ソフトなどを使うことも出来るし、実は「自動微分」という、より高度な方法もあるのだが、まあまずは「既存のライブラリに頼ってみる」能力も必要だろう.
今回は、NLsolve という、非線形方程式の近似根を求めるそのものズバリのライブラリ(Package)を使ってみよう.
まずは NLsolve package を呼び出しておこう
このパッケージをインストールしてない人は、先に Pkg.add("NLsolve")
をしておこう.
using NLsolve
さて,われわれは非線形方程式,例えば $f=0$, を解いてその近似根を求めたい. これを行うそのものズバリの命令が,NLsolve package に含まれている nlsolve 命令だ.
そしてこの命令の使い方だが,(プログラムの都合上)ちょっと変わった仕様なので注意が要る. 具体的には,解きたい方程式 $f(x, param) = 0$ (求めたい根は $x$, param は残りの引数=パラメータ)に対して,
nlsolve( 関数g(配列引数1, 配列引数2), 初期値配列, オプション)
として使うのだが,以下のように注意点がいくつかある.
- 関数 g の引数形式は決まっていて,2つのベクトル値変数を引数に取る.スカラー問題の場合,ベクトルに直す必要がある.
- g の引数の 2つのベクトル値変数を vout, vin という並び1のものとすると,この二つの変数の間には vout $= f$ (vin, param) という関係を設定する必要がある.
- 関数 g 自身の出力はなんでも良い.
- 初期値も必ずベクトル値.
- これの出力はオブジェクトで、見かけは一種の説明文になってしまう.近似根はこの結果の zero 要素で,やはり配列.出力の最後に .zero とつけると取り出せる.
よって,具体的には, 無名関数と名前付き引数の機能を用いて,次のようにするのがまずは簡単なおすすめということになる.
スカラー問題の場合:
解きたい方程式がスカラー変数 $x$ に対してスカラー値関数 $f$ を用いて $f(x, param) = 0$ (param はパラメータ) となっているとして,
nlsolve( (vout,vin)->vout[1]=f(vin[1],param), [ xの初期値 ] ).zero[1]
とすれば近似根を取り出せる2.
ベクトル問題の場合:
解きたい方程式がベクトル値変数 $\boldsymbol{x}$ に対してベクトル値関数 $\boldsymbol{f}$ を用いて $\boldsymbol{f}(\boldsymbol{x}, param) = \boldsymbol{0}$ (param はパラメータ) となっているとして, .= によってベクトルの内容を全コピーできる機能に留意しながら
nlsolve( (vout,vin)->vout.=f(vin,param), ベクトルxの初期値 ).zero
とすれば近似根を取り出せる3.
もっと簡単に ver.1
(ちょっと下に、もうちょっと実用的な ver.2 があるので、実際に使う前にそこもみておくといいんじゃないかな)
上の形でも充分に複雑で面倒なので、もっと簡単にしよう.
それには、using NLsolve
した後に、
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
return nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini]).zero[1]
else
return nlsolve((vout,vin)->vout .= func(vin,params...), ini).zero
end
end
などとお手軽関数 nls
を定義しておいて、これを使えば良い4.
これなら、
スカラー問題の場合: 解きたい方程式がスカラー変数 $x$ に対してスカラー値関数 $f$ を用いて $f(x, params) = 0$ (params はパラメータ, つまり、関数 $f$ の2つ目以降の引数) となっているとして,
nls( f, params, ini = xの初期値 )
とすれば近似根を取り出せる.
文法的に、初期値設定 ini =
は省略できないが、params
はあってもなくても複数でも良い.
ベクトル問題の場合: 解きたい方程式がベクトル値変数 $\boldsymbol{x}$ に対してベクトル値関数 $\boldsymbol{f}$ を用いて $\boldsymbol{f}(\boldsymbol{x}, params) = \boldsymbol{0}$ (params はパラメータ、つまり、関数 $f$ の2つ目以降の引数) となっているならば、
nls( f, params, ini = xの初期ベクトル )
とすれば近似根ベクトルを取り出せる.
文法的に、初期値設定 ini =
は省略できないが、params
はあってもなくても複数でも良い.
スカラー問題でやってみる
まず、パラメータがない、とても簡単な例をやってみよう. $\cos$ 関数の根、つまり、$\cos(x) = 0$ を満たす $x$ を一つ求めてみよう. それには、例えば近似根の初期値を $x_0 = 1.0$ などとして次のようにすれば良い.
nls(cos, ini = 1.0) # ini = 1.0 は近似根 x の探索初期値.
すると
1.570796326795489
とたしかに $\pi/2$ 相当が得られている.
次に、パラメータ(関数 $f$ の2つ目以降の引数のこと)が有る例をやってみよう. 例として、パラメータ $p$ の平方根を求めるような問題を考えてみよう. それには、 関数を $f(x,p) = x^2 - p$ としてこの根 $x$ を求めればいい. だから、次のようになるだろう.
f(x,p) = x^2 - p
nls(f, 2.0, ini=1.0) # 二つ目の引数 2.0 はパラメータ.ini = 1.0 は近似根 x の探索初期値.
すると
1.4142135623746899
と得られて、確かに $\sqrt{2}$ 相当が得られている.
パラメータは複数でも良いので、それを例でみてみよう. 例えば、上の関数を少しだけいじって $f(x,p,q) = x^q - p$ としてこの根 $x$ を求める問題を考えよう.
f(x,p,q) = x^q - p
nls(f, 2.0, 2.1, ini=1.0) # x^(2.1) - 2.0 = 0 の近似根を求める
とすると、
1.391065619249075
を得る. これが確かに近似根になっていることは、
f(1.391065619249075, 2.0, 2.1)
としてみると
9.798828415341632e-12
となるので、確認できる.
もう一つ、似たような関数 $f(u,U)$ に対して $f(u,U)=0$であるような近似根 $u$ を求めてみよう. この場合も $U$ は問題のパラメータ扱いだ.
# 対象の関数を定義.これが ゼロになるような u を求めたい.
Δt, γ = 0.1, 1.0
f(x, U) = x - U - Δt*γ*(sin((x+U)/2))^1.2
nls(f, 3.0, ini = 1.0) # 二つ目の引数 3.0 はパラメータ.ini = 1.0 は近似根 x の探索初期値.
とすると、次のような結果を得る.
3.0091718885079293
さてこの結果が確かに $f \cong 0$ を満たすのか、確認しておこう.
x = nls(f, 3.0, ini=1.0)
f(x, 3.0)
とすると、
5.724587470723463e-17
となり、大変良い精度がでているようだ.うん、悪くない.
ベクトル値問題もやってみよう
では次に,ベクトル問題もやってみよう.
例として,パラメータ $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$ を次のように定義しよう.
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
試しに,適当な値を入れてチェックしてみる.
h([4.0,4.0], 0.5, 0.4)
とすると
2-element Array{Float64,1}:
0.5
-0.2545
となり、まあ良さそうだ.
では,上のご託宣に従って,近似根を求めてみよう.
p = 0.5
q = 0.4
ini_v = [4.0,4.0]
z = nls(h, p, q, ini = ini_v) # パラメータが複数でも大丈夫.
とすると
2-element Array{Float64,1}:
4.48688
3.87691
と得られる. さらに実際にこの値が近似根になっていることも確認しておこう.
h(z, p, q)
とすると
2-element Array{Float64,1}:
1.36335e-13
-1.36668e-13
確かにほぼ ゼロベクトルが出力されるので、OK と言って良さそうだ.
もっと簡単に ver.2
簡単に使えるのはいいけれども、実際にまっとうな近似根が得られているかどうかの「上手く行ったか?という情報」も一緒に欲しいなあということも多いと思う.
自分で確認してもいいけれども、nlsolve
がそれを知っているはずなので、それを教えてもらってもいいだろう.
そういう人は、上の nls
の定義ではなくて、ちょっと変えた、たとえば下のような定義を用いるのがいいんじゃないかな.
これは、using NLsolve
した後に、
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
v = r.zero[1]
else
r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
v = r.zero
end
return v, r.f_converged
end
とすればいい5.
これだと、近似根と、「その近似根で $f \cong 0$ になっているか?」という二つの情報が出力される.
ちなみに、$f \cong 0$ とみなす判定基準は nlsolve
のオプション ftol =
で変えられる(ftol
のデフォルトは 1.0e-8 だね).
試しに、上の $\cos$ の例で
nls(cos, ini = 1.0)
としてみると、
(1.570796326795489,true)
と出力されて、一つ目の数字が近似根で、2つめの “true” が「この近似根で $f \cong 0$ となっている」ことがわかる.
実際、根が求まらない単純なケースとして、負の数字の平方根を求めることに相当する近似根計算をさせてみよう.
f(x,p) = x^2 - p
nls(f, -0.0001, ini = 1.0)
とすると、
(-4.385150012653803e-11,false)
と出力され、(探索ルーチンが頑張って求めた)一見近似根っぽい数字が一つ目の出力項に出ているけれども、2つめの出力項の “false” によって、それは $f \cong 0$ を満足するような近似根ではない、という情報がわかる.
「ややこしい関数で近似根を求めさせると失敗するかも」というようなケースでは、こうした情報をチェックするようなプログラムを書いたほうが良さそうだね.