Octave で数値解析のために必要なプログラミングと，いくつかの応用例を示します。
Octave の文末は，単に改行するか，; で終わって改行するかです。; で終わる場合は結果が表示されません。
1 + 2
1 + 2;
条件分岐の if 文の書式は以下のようになっています。
if (条件式1)
条件式1 が満たされた場合に実行する文
elseif (条件式2)
条件式2 が満たされた場合に実行する文
else
それ以外の場合に実行する文
endif
% rand() は (0, 1) 区間のランダムな実数を生成する。
x = rand();
if (x <= 0.5)
disp(x), disp(" は 0.5 以下")
else
disp(x), disp(" は 0.5 より大きい")
endif
$\displaystyle \sum_{i=1}^{5} i\ $ を計算する例
isum = 0;
for i = 1:5
isum = isum + i;
endfor
disp(isum)
isum = 0;
i = 1;
while (i <= 5)
isum = isum + i;
i = i + 1;
endwhile
disp(isum)
$\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$
function sumi = isum(n)
if (n == 1)
sumi = 1;
else
sumi = isum(n-1) + n;
endif
endfunction
disp(isum(5))
sum()¶Octave にはベクトル要素の和を求める関数 sum() が組み込まれています。
sum([1:5])
解析的な微分の定義は（前方差分の極限として）
$$\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 を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。
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 では必要でした。）
%plot --format svg
x = [-2*pi:4*pi/100:2*pi];
plot(x, cos(x) - diff_f(x, 1.0E-6))
以下の例からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけでもなさそうです。（縦軸の目盛りに注目。）
また，デフォルトでは linewidth が 0.5 らしく，細すぎに見えるので，1 に設定してみます。
plot(x, cos(x) - diff_f(x, 1.0E-7), "linewidth", 1)
シンプソン法で $\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 文の繰り返し処理で和をとっています。
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)
quad() 等の利用¶Octave には数値積分をしてくれる quad() 等の関数が用意されています。用意されている関数を使うと楽です。
format long g;
quad("f", 0, pi/2)
quadl("f", 0, pi/2)
quadv("f", 0, pi/2)
quadcc("f", 0, pi/2)
$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。
function y = sinsin(x)
y = sin(sin(x));
endfunction
quad("sinsin", 0, pi)
quadl("sinsin", 0, pi)
$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$ の範囲でグラフを描きます。
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;
while 文による解法¶上のグラフから，$f(x) = 0$ となる $x$ の1つは，$2 < x < 3$ にあることがわかります。
探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件は，数値解の精度に関係しますので，いくつか値を変えて出力させます。
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 を変えても変わらない桁までが，精度のよい数値解ということになります。上の例では，
2.404825557695586
くらいでしょうか。
f(2.404825557695586)
fsolve() の利用¶Octave には非線形方程式を解いてくれる関数 fsolve() が用意されています。
$f(x) = 0$ を $x_0 = 2.5$ を探索の初期値として数値解を求める例が以下です。
fsolve("f", 2.5)
f(ans)
fsolve() はデフォルトでは 1e-6 程度の精度で求めるようです。許容誤差をもっと小さくして求めてみます。
fsolve("f", 2.5, optimset("TolX", 1.E-14, "TolFun", 1.E-14))
f(ans)
4次のルンゲ・クッタ法で常微分方程式を解きます。
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$ 分割して解きます。
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 を変えて結果をみます。分割数を変えても変わらない桁数までをもって，精度の良い数値解とします。
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 を変えても変わらないところまでというと，
0.99959156
くらいが最後の値になるでしょうか。
plot¶plot(T, X, "linewidth", 1)
xlim([0 10]); ylim([0 1]);
xlabel("t"); ylabel("x");
grid on;
lsode() の利用¶Octave には常微分方程式を数値的に解いてくれる関数 lsode() が用意されています。それを使ってみましょう。（すでに定義されている函数を使うほうが楽ちん。）
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)
デフォルトだと少し精度が足りないかも知れません。
lsode() の精度がどうなっているか確かめます。
lsode_options()
absolute tolerance と relative tolerance の値を以下のように設定してもう一度計算してみます。
lsode_options("absolute tolerance", 1e-14)
lsode_options("relative tolerance", 1e-14)
t = linspace(t0, t1, N);
x = lsode("f", x0, t);
x(N)
次のような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) $$を，
して解きます。
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)]$ というイメージです。）
function dx = F(x, t, p)
dx(1) = x(2);
dx(2) = -x(1) - p(1)*x(2) + p(2)*cos(t);
endfunction
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 の渡し方
plot(t, x0000(:, 1), "linewidth", 1,
t, x0502(:, 1), "linewidth", 2)
legend("a = 0, b = 0", "a = 0.5, b = 0.2");
grid("on");
上の計算例では，lsode() 関数を使って解いた数値解は x0000 や x0502 に N行2列の行列の形で格納されています。
size(x0502)
1行目を取り出して plot するときに x0000(:, 1) と書いています。
行列の要素（特定の列とか特定の行とか特定の成分）を取り出す例を以下に示します。
M は3行2列の行列です。
M = [1 2; 3 4; 5 6]
size(M)
M の1列目を取り出すには以下のようにします。
M(:, 1)
M の2行目を取り出すには...
M(2, :)
M の1行2列目の成分を取り出すには...
M(1, 2)