はじめに:Octave で数値解析

Octave で数値解析のために必要なプログラミングと,いくつかの応用例を示します。

Octave の文末は,単に改行するか,; で終わって改行するかです。; で終わる場合は結果が表示されません。

In [1]:
1 + 2
ans = 3
In [2]:
1 + 2;

条件分岐

条件分岐の if 文の書式は以下のようになっています。

if (条件式1)
    条件式1 が満たされた場合に実行する文
elseif (条件式2)
    条件式2 が満たされた場合に実行する文
else
    それ以外の場合に実行する文
endif
In [3]:
% rand() は (0, 1) 区間のランダムな実数を生成する。
x = rand();
if (x <= 0.5)
    disp(x), disp(" は 0.5 以下")
else
    disp(x), disp(" は 0.5 より大きい")
endif
0.047174
 は 0.5 以下

繰り返し処理

$\displaystyle \sum_{i=1}^{5} i\ $ を計算する例

for

for 文の書式は以下の通りです。

for i = 開始値:増分:終端値
    i が開始値以上,終端値以下の場合に実行される文
endfor
In [4]:
isum = 0;
for i = 1:5
    isum = isum + i;
endfor
disp(isum)
15

while

while 文の書式は以下の通りです。

while (条件式)
    条件式が成り立つ場合に実行する式
endwhile
In [5]:
isum = 0;
i = 1;
while (i <= 5)
    isum = isum + i;
    i = i + 1;
endwhile
disp(isum)
15

再帰的定義関数

$\displaystyle \mbox{isum}(n) = \sum_{i=1}^{n} i , \ (n \geq 1)\ $ を以下のように考えます。

もし,$n = 1$ なら $\mbox{isum}(n) = 1$

もし,$n > 1$ なら $\mbox{isum}(n) = \mbox{isum}(n-1) + n$

In [6]:
function sumi = isum(n)
    if (n == 1)
        sumi = 1;
    else 
        sumi = isum(n-1) + n;
    endif
endfunction
In [7]:
disp(isum(5))
15

参考:sum()

Octave にはベクトル要素の和を求める関数 sum() が組み込まれています。

In [8]:
sum([1:5])
ans = 15

数値微分

解析的な微分の定義は(前方差分の極限として)

$$\frac{df}{dx} \equiv \lim_{h \rightarrow 0} \frac{f(x+h) - f(x)}{h}$$

でした。滑らかな関数であれば,(後方差分の極限として)

$$\frac{df}{dx} = \lim_{h \rightarrow 0} \frac{f(x) - f(x-h)}{h}$$

としても良いです。

現実世界では,$ h\rightarrow 0$ の極限はとれませんから十分小さい値として h を定義し,近似的な数値微分を(前方差分と後方差分の平均である中心差分として)以下のように定義しましょう。

\begin{eqnarray} \frac{df}{dx} &\simeq& \frac{1}{2} \left\{\frac{f(x+h) - f(x)}{h} + \frac{f(x) - f(x-h)}{h} \right\} \\ &=& \frac{f(x+h) - f(x-h)}{2h} \end{eqnarray}
In [9]:
function ret = f(x)
    ret = sin(x);
endfunction

function diff = diff_f(x, h)
    diff = (f(x + h) - f(x - h))/(2*h);
endfunction

以下の例では,$f(x) = \sin x$ に対して $f'(x) = \cos x$ と,上記の数値微分との誤差を描いています。

また,必要なら以下のようにフォーマットを SVG に設定します。 (brew install octave した macOS では必要でした。)

In [12]:
%plot --format svg
In [13]:
x = [-2*pi:4*pi/100:2*pi];
plot(x, cos(x) - diff_f(x, 1.0E-6))
Gnuplot Produced by GNUPLOT 5.4 patchlevel 1 -2e-10 -1e-10 0 1e-10 2e-10 -10 -5 0 5 10 gnuplot_plot_1a

以下の例からわかるように,h をやみくもに小さくすればするほど,数値微分の精度があがる,というわけでもなさそうです。(縦軸の目盛りに注目。)

また,デフォルトでは linewidth が 0.5 らしく,細すぎに見えるので,1 に設定してみます。

In [14]:
plot(x, cos(x) - diff_f(x, 1.0E-7), "linewidth", 1)
Gnuplot Produced by GNUPLOT 5.4 patchlevel 1 -4e-09 -3e-09 -2e-09 -1e-09 0 1e-09 2e-09 3e-09 -10 -5 0 5 10 gnuplot_plot_1a

数値積分

シンプソン法で $\displaystyle \int_a^b f(x)\, dx$ を求める。

積分区間 $[a, b]$ を $N (=2n)$ 等分(偶数等分)し,

$$h = \frac{b-a}{N}, \quad f_i = f(a + i\,h)$$

とする。シンプソン法の公式は

\begin{eqnarray} \int_a^b f(x)\, dx &\simeq& \frac{h}{3}( f_0 + 4 f_1 + f_2) + \frac{h}{3}( f_2 + 4 f_3 + f_4) + \cdots\\ &&\quad \cdots + \frac{h}{3}( f_{N-2} + 4 f_{N-1} + f_N)\\ &=&\frac{h}{3} \left(f_0 + f_{2n} + 2 \sum_{i=1}^{n-1} f_{2i} + 4 \sum_{i=1}^{n} f_{2n-1}\right) \end{eqnarray}

for 文による解法

$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を,積分区間を $N$ 分割してシンプソン法で解く例です。精度確認のため,分割数を2倍にしてもう一度計算しています。

for 文の繰り返し処理で和をとっています。

In [15]:
format long g;
function y = f(x)
    y = sin(x);
endfunction

a = 0;    % 積分の下端
b = pi/2; % 積分の上端
N = 20;   % 分割数 N は偶数とすること(シンプソン法の決め事)
h = (b - a)/N;

simpson_f = 0;
for i = 2:2:N
    simpson_f  = (simpson_f  + 
                  h/3*(f(a+(i-2)*h) + 4*f(a+(i-1)*h) + f(a+i*h)));
endfor
printf("%d 分割の数値解は ", N), disp(simpson_f)

% 分割数を倍にしてもう一度計算
N = 40;
h = (b - a)/N;

simpson_f = 0;
for i = 2:2:N
    simpson_f  = (simpson_f  + 
                  h/3*(f(a+(i-2)*h) + 4*f(a+(i-1)*h) + f(a+i*h)));
endfor
printf("%d 分割の数値解は ", N), disp(simpson_f)
20 分割の数値解は 1.000000211546591
40 分割の数値解は 1.000000013214379

参考:quad() 等の利用

Octave には数値積分をしてくれる quad() 等の関数が用意されています。用意されている関数を使うと楽です。

In [16]:
format long g;
quad("f", 0, pi/2)
quadl("f", 0, pi/2)
quadv("f", 0, pi/2)
quadcc("f", 0, pi/2)
ans = 0.9999999999999999
ans =                1
ans = 0.9999999978731193
ans = 0.9999999999999998

$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。

In [17]:
function y = sinsin(x)
    y = sin(sin(x));
endfunction

quad("sinsin", 0, pi)
quadl("sinsin", 0, pi)
ans = 1.786487481950052
ans = 1.786487481953686

方程式の数値解

$f(x) = 0$ の解を求める。

まず,何らかの方法で(例えば,plot f(x) としてグラフを描いてみて $x$ 軸と交差する点のあたりをつけるなどして)$f(x) = 0$ の解である $x$ に近いと思われる値 $x_k$ がわかったとする。

この $x_k$ のまわりに $f(x)$ をテイラー展開すると,

$$ f(x) \simeq f(x_k) + f'(x_k)\,(x - x_k)$$

$f(x) = 0$ であるから,上式の左辺をゼロとおいて,$x$ について解くと,

$$ x \equiv x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$

これを $|f(x_{k+1})| < \epsilon$ (ここで $\epsilon$ 十分小さい数値,たとえば $1.0\times 10^{-8}$ とか?)になるまで反復計算を行い,このときの $x_{k+1}$ をもって,$f(x) = 0$ の近似解とする。

$f'(x_k)$ を求めるのが面倒な場合は,以下のようにして数値的な微分で代用する。

$$f'(x_k) \simeq \frac{f(x_k + h) - f(x_k - h)}{2 h}, \quad |h| \ll 1$$

つまり,

$$x_{k+1} = x_k - \frac{f(x_k)\times 2 h}{f(x_k + h) - f(x_k - h)}$$

例:ベッセル関数のゼロ点

以下では例としてベッセル関数 besselj(0, x) のゼロ点を求めています。

まず,$0 \leq x \leq 10$ の範囲でグラフを描きます。

In [18]:
function y = f(x)
    y = besselj(0, x);
endfunction

x = [0:0.1:10];
plot(x, f(x), "linewidth", 1)
legend("J_0 Bessel funcion");
grid on;
Gnuplot Produced by GNUPLOT 5.4 patchlevel 1 -1 -0.5 0 0.5 1 0 2 4 6 8 10 J0 Bessel funcion J0 Bessel funcion

while 文による解法

上のグラフから,$f(x) = 0$ となる $x$ の1つは,$2 < x < 3$ にあることがわかります。

探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件は,数値解の精度に関係しますので,いくつか値を変えて出力させます。

In [19]:
x = 2.5;

h = 1.E-8
while (abs(f(x)) > h)
       x = x - f(x)*2*h/(f(x + h)-f(x - h));
endwhile
disp(x)

h = 1.E-10
while (abs(f(x)) > h)
       x = x - f(x)*2*h/(f(x + h)-f(x - h));
endwhile
disp(x)

h = 1.E-12
while (abs(f(x)) > h)
       x = x - f(x)*2*h/(f(x + h)-f(x - h));
endwhile
disp(x)
h =            1e-08
2.404825557695586
h =            1e-10
2.404825557695586
h =            1e-12
2.404825557695586

解の収束条件 h を変えても変わらない桁までが,精度のよい数値解ということになります。上の例では,

2.404825557695586

くらいでしょうか。

In [20]:
f(2.404825557695586)
ans = 9.701672147389553e-14

参考:fsolve() の利用

Octave には非線形方程式を解いてくれる関数 fsolve() が用意されています。

$f(x) = 0$ を $x_0 = 2.5$ を探索の初期値として数値解を求める例が以下です。

In [21]:
fsolve("f", 2.5)
f(ans)
ans = 2.404824591812223
ans = 5.014361283104874e-07

fsolve() はデフォルトでは 1e-6 程度の精度で求めるようです。許容誤差をもっと小さくして求めてみます。

In [22]:
fsolve("f", 2.5, optimset("TolX", 1.E-14, "TolFun", 1.E-14))
f(ans)
ans = 2.404825557695772
ans = 1.535884772676773e-16

微分方程式の数値解法

4次のルンゲ・クッタ法で常微分方程式を解きます。

1階常微分方程式

1階の常微分方程式

$$ \frac{dx}{dt} = f(x, t)$$

を4次のルンゲ・クッタ法で解く例です。

初期条件を $t=t_0$ で $x(t_0) = x_0$ とし,$t = t_1$ までを $N$ 分割して解きます。

$$ h = \frac{t_1 - t_0}{N}, \quad t_i = t_0 + i\, h, \quad x_i = x(t_i)$$

とし,以下の式から次のステップの値 $x_{i+1}$ を決めます。

\begin{eqnarray} k_1 &=& h\, f(x_i, t_i) \\ k_2 &=& h\, f(x_i + k_1/2, t_i + h/2) \\ k_3 &=& h\, f(x_i + k_2/2, t_i + h/2) \\ k_4 &=& h\, f(x_i + k_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \end{eqnarray}

例:ロジスティック方程式

例として,簡単なロジスティック方程式

$$ \frac{dx}{dt} = (1-x)\,x$$

を,初期条件を $t_0 = 0$ で $x_0 = x(t_0) = 0.1$ とし,$t = t_0 = 0$ から $t_1 = 10$ までを $ N = 100$ 分割して解きます。

In [23]:
function y = f(x, t)
    y = (1 - x) * x;
endfunction

t0 = 0;
x0 = 0.1;
t1 = 10;
N = 100;

h = (t1 - t0)/N;

t = t0;
x = x0;
T = [];
X = [];
T(1) = t;
X(1) = x;

for i = 2:N+1
    k1 = h*f(x, t);
    k2 = h*f(x+0.5*k1, t+0.5*h);
    k3 = h*f(x+0.5*k2, t+0.5*h);
    k4 = h*f(x+k3, t+h);
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6;
    t = t + h;
    T(i) = t;
    X(i) = x;
endfor

printf("N = %d   X(N+1) = %18.15f \n", N, X(N+1))
N = 100   X(N+1) =  0.999591565221868 

分割数 N を変えて結果をみます。分割数を変えても変わらない桁数までをもって,精度の良い数値解とします。

In [24]:
printf("N = %d   X(N+1) = %18.15f \n", N, X(N+1))

t0 = 0;
x0 = 0.1;
t1 = 10;
N = 400;

h = (t1 - t0)/N;

t = t0;
x = x0;
T = [];
X = [];
T(1) = t;
X(1) = x;

for i = 2:N+1
    k1 = h*f(x, t);
    k2 = h*f(x+0.5*k1, t+0.5*h);
    k3 = h*f(x+0.5*k2, t+0.5*h);
    k4 = h*f(x+k3, t+h);
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6;
    t = t + h;
    T(i) = t;
    X(i) = x;
endfor

printf("N = %d   X(N+1) = %18.15f \n", N, X(N+1))
N = 100   X(N+1) =  0.999591565221868 
N = 400   X(N+1) =  0.999591567508863 

N を変えても変わらないところまでというと,

0.99959156

くらいが最後の値になるでしょうか。

データの plot

In [25]:
plot(T, X, "linewidth", 1)
xlim([0 10]); ylim([0 1]);
xlabel("t"); ylabel("x");
grid on;
Gnuplot Produced by GNUPLOT 5.4 patchlevel 1 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 x t gnuplot_plot_1a

参考:lsode() の利用

Octave には常微分方程式を数値的に解いてくれる関数 lsode() が用意されています。それを使ってみましょう。(すでに定義されている函数を使うほうが楽ちん。)

In [26]:
function y = f(x, t)
    y = (1 - x) * x;
endfunction

t0 = 0;
x0 = 0.1;
t1 = 10;
N = 100;

t = linspace(t0, t1, N);
x = lsode("f", x0, t);
format long
x(N)
ans = 0.999591534534071

デフォルトだと少し精度が足りないかも知れません。

lsode() の精度がどうなっているか確かめます。

In [27]:
lsode_options()
Options for LSODE include:

  keyword                                             value
  -------                                             -----
  absolute tolerance                                  1.49012e-08
  relative tolerance                                  1.49012e-08
  integration method                                  stiff
  initial step size                                   -1
  maximum order                                       -1
  maximum step size                                   -1
  minimum step size                                   0
  step limit                                          100000

absolute tolerancerelative tolerance の値を以下のように設定してもう一度計算してみます。

In [28]:
lsode_options("absolute tolerance", 1e-14)
lsode_options("relative tolerance", 1e-14)
In [29]:
t = linspace(t0, t1, N);
x = lsode("f", x0, t);

x(N)
ans = 0.999591567517052

2階常微分方程式

次のような2階常微分方程式を解く例を示します。

$$ \frac{d^2 x}{dt^2 } = f\left(x, \frac{dx}{dt}, t\right)$$

この場合には, $\displaystyle x_1 \equiv x, x_2 \equiv \frac{dx}{dt}$ とおけば,次のような連立1階常微分方程式の形に帰着できます。

\begin{eqnarray} \frac{dx_1}{dt} &=& F_1(x_1, x_2, t) = x_2 \\ \frac{dx_2}{dt} &=& F_2(x_1, x_2, t) = f(x_1, x_2, t) \end{eqnarray}

連立であれ何であれ,1階常微分方程式の形になったら,ルンゲ・クッタ法で解くことができます... が,すでに Octave の lsode() 関数を使っているので,これを利用して連立1階常微分方程式を解いてみます。

例:減衰+強制振動

例として,簡単な減衰+強制振動の方程式

$$ \frac{d^2 x}{dt^2} = -x - a \frac{dx}{dt} + b \cos(t) $$

を,

  • $a , b $ にいくつか値を入れて
  • 初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 3, v_0 = v(t_0) = 0$,
  • $t = t_0 = 0$ から $t_1 = 20$ までを $ N $ 分割

して解きます。

lsode() で連立1階方程式

まず,$\displaystyle x_1 \equiv x, x_2 \equiv \frac{dx}{dt}$ として,連立1階常微分方程式の形にします。

また,変数 $x_1, x_2$ 以外にも $a, b$ をパラメータとして含むので,$p_1 \equiv a, p_2 \equiv b$ とします。

\begin{eqnarray} \frac{dx_1}{dt} &=& F_1(x_1, x_2, t, p_1, p_2) = x_2 \\ \frac{dx_2}{dt} &=& F_2(x_1, x_2, t, p_1, p_2) = -x_1 - p_1 x_2 + p_2 \cos t \end{eqnarray}

Octave の lsode() 関数を使うためには以下のように定義します。

($\boldsymbol{x} = [x(1), x(2)], \boldsymbol{F} = [F(1), F(2)], \boldsymbol{p} = [p(1), p(2)]$ というイメージです。)

In [30]:
function dx = F(x, t, p)
    dx(1) = x(2);
    dx(2) = -x(1) - p(1)*x(2) + p(2)*cos(t);
endfunction
In [31]:
t0 = 0;
x0 = [3 0];  % 初期条件:x(1) = x = 3, x(2) = v = 0
t1 = 20;
N = 100;

P = [0 0];       % a = 0, b = 0
t = linspace(t0, t1, N);
x0000 = lsode(@(x, t)F(x, t, P), x0, t); % パラメータ P の渡し方

P = [0.5 0.2];       % a = 0.5, b = 0.2
t = linspace(t0, t1, N);
x0502 = lsode(@(x, t)F(x, t, P), x0, t); % パラメータ P の渡し方

2つのグラフを重ねて表示

In [32]:
plot(t, x0000(:, 1), "linewidth", 1, 
     t, x0502(:, 1), "linewidth", 2)
legend("a = 0, b = 0", "a = 0.5, b = 0.2");
grid("on");
Gnuplot Produced by GNUPLOT 5.4 patchlevel 1 -3 -2 -1 0 1 2 3 0 5 10 15 20 a = 0, b = 0 a = 0, b = 0 a = 0.5, b = 0.2 a = 0.5, b = 0.2

参考:行列の要素を取り出す

上の計算例では,lsode() 関数を使って解いた数値解は x0000x0502 に N行2列の行列の形で格納されています。

In [33]:
size(x0502)
ans =

   100     2

1行目を取り出して plot するときに x0000(:, 1) と書いています。

行列の要素(特定の列とか特定の行とか特定の成分)を取り出す例を以下に示します。

M は3行2列の行列です。

In [34]:
M = [1 2; 3 4; 5 6]
size(M)
M =

   1   2
   3   4
   5   6

ans =

   3   2

M の1列目を取り出すには以下のようにします。

In [35]:
M(:, 1)
ans =

   1
   3
   5

M の2行目を取り出すには...

In [36]:
M(2, :)
ans =

   3   4

M の1行2列目の成分を取り出すには...

In [37]:
M(1, 2)
ans = 2