葛西 真寿 弘前大学大学院理工学研究科
Maxima はコンピュータ代数システム（いわゆる数式処理システム）ですが，構造化プログラミングの機能もあり，条件分岐や繰り返し処理などを利用することで，プログラミング言語としても利用できます。ここでは，プログラミングに必要なコマンドと，いくつかの数値解析への応用例をまとめておきます。
Maxima の文末は，; か $ で終わります。$ で終わる場合は結果が表示されません。
1 + 2;
1 + 2$
書式
if 条件式1 then
条件式1が満たされる場合に実行する文
elseif 条件式2 then
条件式2が満たされる場合に実行する文
else
上記以外の場合に実行する文 $x: random(1.0)$
if x < 0.5 then
print(x, " は 0.5 より小さい")
elseif x > 0.7 then
print(x, " は 0.7 より大きい")
else
print(x, " は 0.5 以上 0.7 以下")$
$\displaystyle \sum_{i = 1}^{5} i$ を計算する例
for ... do¶書式
for i: istart thru iend step istep do(実行文1, 実行文2, ...)sum: 0$
for i: 1 thru 5
do(sum: sum + i,
print("1 から ", i, " までの和は ", sum)
)$
while ... do¶sum: 0$
i: 1$
while i <= 5
do(sum: sum + i,
print("1 から ", i, " までの和は ", sum),
i: i + 1
)$
$\mbox{isum}(n) = \displaystyle \sum_{i=1}^{n} i$ を以下のように考えます。
もし，$n = 1$ なら $\mbox{isum}(n) = \mbox{isum}(1) = 1$
もし，$n > 1$ なら $\mbox{isum}(n) = \mbox{isum}(n-1) + n$
isum(n):= if n = 1 then 1 else isum(n-1) + n $
/* i = 1 から 5 まで増分 2 ごとに表示する例 */
for i:1 thru 5 step 2
do(print("1 から ", i, " までの和は ", isum(i)))$
sum, nusum¶上の例で，isum(n) という関数を再帰的に定義しましたが，定義などしなくても，実は数列の和を計算する関数が Maxima には組み込まれています。
sum(i, i, 1, 5);
Maxima はコンピュータ代数システムですから，一般の $n$ についても数列の和を計算してくれそうですが，sum() だと...
sum(k, k, 1, n);
組込関数 nusum() なら，数列の和の公式を知っています。
nusum(k, k, 1, n);
sum(k**2, k, 1, n) = nusum(k**2, k, 1, n);
解析的な微分の定義は（前方差分の極限として） $$\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}$$ としても良いです。
現実世界の gnuplot では，$ h\rightarrow 0$ の極限はとれませんから十分小さい値として h を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。
以下の例では，$f(x) = \sin x$ に対して $f'(x) = \cos x$ と，上記の数値微分との誤差を描いています。
/* Jupyter Notebook にインラインでグラフを表示させる場合 */
set_plot_option([svg_file, "~/.maxplot.svg"])$
kill(x)$
f(x):= sin(x)$
/* f(x) の数値微分 */
diff_f(x, h):= (f(x+h) - f(x-h))/(2.0*h)$
plot2d(cos(x) - diff_f(x, 1.0E-6), [x, -2*%pi, 2*%pi], [ylabel,""])$
以下の例からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけでもなさそうです。
plot2d(cos(x) - diff_f(x, 1.0E-7), [x, -2*%pi, 2*%pi], [ylabel,""])$
diff¶実は数値微分などしなくても，Maxima は関数の微分をしてくれます。
f(x):= sin(x)$
diff(f(x), x);
シンプソン法で $\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 ... do 文による解法¶$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を，積分区間を $N0$ 分割してシンプソン法で解く例です。
for ... do 文の繰り返し処理で和をとっています。
f(x):= sin(x)$ /* 被積分関数 */
a: 0$
b: %pi/2$
N: 10$
h: float((b-a)/N)$
xi(i):= a + i*h$
fi(i):= f(xi(i))$
simpson_f: 0$
for i: 2 thru N step 2
do(simpson_f: simpson_f + h/3*(fi(i-2) + 4*fi(i-1) + fi(i)))$
print(N, " 分割の数値解は ", simpson_f)$
/* 精度確認のため，分割数を倍にして計算 */
N: 2*N$
h: float((b-a)/N)$
xi(i):= a + i*h$
fi(i):= f(xi(i))$
simpson_f: 0$
for i: 2 thru N step 2
do(simpson_f: simpson_f + h/3*(fi(i-2) + 4*fi(i-1) + fi(i)))$
print(N, " 分割の数値解は ", simpson_f)$
kill(a, b, N, i, h)$
f(x):= sin(x)$
/* a 積分の下端 */
/* b 積分の上端 */
/* N 分割数 N は偶数とすること（シンプソン法の決め事） */
simp_f(a, b, N, i):=
block([h, xi],
h: float((b-a)/N), xi(i):= float(a + i*h),
if i = 2 then
h/3*(f(xi(i-2)) +4*f(xi(i-1)) + f(xi(i)))
else
simp_f(a, b, N, i-2) +
h/3*(f(xi(i-2)) +4*f(xi(i-1)) + f(xi(i)))
)$
/* 実際の計算の際は，simp_f(a, b, N, N) と同じ N を2回書く */
N: 10$
print(N, " 分割の数値解は ", simp_f(0, %pi/2, N, N))$
N: 20$
print(N, " 分割の数値解は ", simp_f(0, %pi/2, N, N))$
integrate, romberg¶実は数値積分などしなくても，Maxima は解析的に積分してくれます。
integrate(sin(x), x, 0, %pi/2);
$\displaystyle \int_0^{\pi} \sin (\sin x) \, dx$ のように解析的に積分できない場合は，romberg という組込関数で数値積分できます。
integrate(sin(sin(x)), x, 0, %pi);
romberg(sin(sin(x)), x, 0, %pi);
数値積分 romberg の精度は，大域変数 rombergtol で決まります。
rombergtol の値を小さくすることで精度をあげることができます。
デフォルトの rombergtol の値は...
rombergtol;
rombergtol の値を変えて計算し，精度を確認します。
rombergtol: 1.e-6$
romberg(sin(sin(x)), x, 0, %pi);
rombergtol: 1.e-8$
romberg(sin(sin(x)), x, 0, %pi);
rombergtol: 1.e-10$
romberg(sin(sin(x)), x, 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)}$$
まず，解を求める関数 $f(x)$ を定義し，$f(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。
以下では例として組み込み関数の1つである bessel_j(0,x) のゼロ点を求めています。Maxima の組み込み関数については，たとえば以下を参照。
f(x):= bessel_j(0,x)$
plot2d(f(x), [x, 0, 10], [legend, "J_0 ベッセル関数"], grid2d)$
while 文による解法¶$f(x) = 0$ となる $x$ の1つは，$2$ から $3$ の間にあることがわかります。
探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。
f(x):= bessel_j(0,x)$
/* x には探索の初期値を入れる */
/* h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。*/
x: 2.5$
h: 1.E-6$
while abs(f(x)) > h
do(x: x - f(x)*2*h/(f(x + h) - f(x - h)))$
x;
/* 精度確認のため，h を変えて計算 */
x: 2.5$
h: 1.E-8$
while abs(f(x)) > h
do(x: x - f(x)*2*h/(f(x + h) - f(x - h)))$
x;
以下のように，再帰定義により $f(x) = 0$ の解を求める関数 newton_f(x, h) を定義します。
探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。
/* 再帰的定義関数による解 */
/* x には探索の初期値を入れる */
/* h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。*/
newton_f(x, h):= if abs(f(x)) > h
then
newton_f(x-f(x)*2*h/(f(x+h)-f(x-h)), h)
else
x $
/* 解の表示 */
newton_f(2.5, 1E-6);
newton_f(2.5, 1E-8);
newton_f(2.5, 1E-10);
newton_f(2.5, 1E-12);
newton_f(2.5, 1E-14);
解の収束条件 h を変えても変わらない桁までが，精度のよい数値解ということになります。上の例では，
2.4048255576957
くらいでしょうか。
f(2.4048255576957);
find_root¶実は，方程式の数値解を求める関数 find_root が Maxima には組み込まれています。
以下は，$2 < x < 3$ の範囲で $f(x) = 0$ の解を数値的に求める例です。
kill(x)$ /* 変数 x に値が入っていると困るので初期化します */
find_root(f(x) = 0, x, 2, 3);
f(%);
上の例で，% は直前の出力（計算結果）を表します。
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$ 分割して解きます。
f(x, t):= (1-x)*x$
t0: 0 $
x0: 0.1 $
t1: 10 $
N: 100 $
h: float((t1 - t0) * 1.0/N)$
t: t0 $
x: x0 $
T[0]: t0$
X[0]: x0$
for i: 1 thru N do(
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[i]: x + (k1 + 2.*k2 + 2.*k3 + k4)/6,
T[i]: t + h,
x: X[i],
t: T[i]
)$
plot2d¶data: makelist([T[i], X[i]], i, 0, N)$
plot2d([discrete, data], [y, 0, 1], grid2d)$
分割数 N を変えて結果をみます。分割数を変えても変わらない桁数までをもって，精度の良い数値解とします。
f(x, t):= (1-x)*x$
t0: 0 $
x0: 0.1 $
t1: 10 $
N: 100$
h: float((t1 - t0) * 1.0/N)$
t: t0 $
x: x0 $
T[0]: t0$
X[0]: x0$
for i: 1 thru N do(
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[i]: x + (k1 + 2.*k2 + 2.*k3 + k4)/6,
T[i]: t + h,
x: X[i],
t: T[i]
)$
print(N, " 分割時の数値 ",X[N])$
N: 400$
h: float((t1 - t0) * 1.0/N)$
t: t0 $
x: x0 $
T[0]: t0$
X[0]: x0$
for i: 1 thru N do(
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[i]: x + (k1 + 2.*k2 + 2.*k3 + k4)/6,
T[i]: t + h,
x: X[i],
t: T[i]
)$
print(N, " 分割時の数値 ",X[N])$
rk¶実は，Maxima にはルンゲ・クッタ法によって常微分方程式を解く関数 rk があらかじめ組み込まれています。使い方は以下のとおりです。
kill(x, t)$
f(x, t):= (1-x)*x$
t0: 0 $
x0: 0.1 $
t1: 10 $
N: 100$
h: float((t1 - t0) * 1.0/N)$
results: rk(f(x, t), x, x0, [t, t0, t1, h])$
plot2d([discrete, results],[y, 0, 1], grid2d)$
関数 rk の計算結果は，上の例ではリスト results に代入されます。
results の行数は 101，1行目が初期条件の [t0, x0]，最終行が [t1, x(t1)] です。
length(results);
results[1];
results[101];
数値解析の本来の目的は，解析的に解けない問題を何とかして解きたい，解いたときの精度もちゃんと把握したい，ということです。
なので，上記の数値解がどの程度の精度で精確であるかを調べるために解析解と比較するのは本来は邪道なのですが，たまたま，上記の微分方程式は解析的に解けて，以下のようになります。
$$x(t) = \frac{\exp(t)}{9 + \exp(t)}$$実は，Maxima には微分方程式を解析的に解く組込関数 ode2 があります。ただし，今回のような場合（変数分離形で解けるタイプ）には，以下のように計算結果がイマイチです。
ode2('diff(x, t) = (1-x)*x, x, t);
ここからは，人間の手で解きましょう。上記の $\%c$ は積分定数なので以下では $C$
まず，左辺は $$ \log |x| - \log|1-x| = \log \left|\frac{x}{1-x}\right|$$ となるので， $$\frac{x}{1-x} = K e^t$$ となることはわかるであろう。ここで，$K \equiv e^C$
初期条件をいれると $$\frac{1/10}{1-1/10} = \frac{1}{9} = K$$ 最終的に $$ x = \frac{K e^t}{1 + K e^t} = \frac{e^t}{9 + e^t}$$となる。
次のような2階常微分方程式をルンゲ・クッタ法で解く。 $$ \frac{d^2 x}{dt^2 } = f\left(x, \frac{dx}{dt}, t\right)$$
この場合には， $\displaystyle v \equiv \frac{dx}{dt}$ とおけば，次のような連立1階常微分方程式の形に帰着できる。 \begin{eqnarray} \frac{dx}{dt} &=& F_1(x, v, t) = v \\ \frac{dv}{dt} &=& F_2(x, v, t) = f(x, v, t) \end{eqnarray}
初期条件を $t=t_0$ で $x(t_0) = x_0, v(t_0) = v_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), \quad v_i = v(t_i)$$とし，以下の式から次のステップの値 $x_{i+1}, v_{i+1}$ を決めます。
\begin{eqnarray} k_1 &=& h\, F_1(x_i, v_i, t_i) \\ m_1 &=& h\, F_2(x_i, v_i, t_i) \\ k_2 &=& h\, F_1(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ m_2 &=& h\, F_2(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ k_3 &=& h\, F_1(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ m_3 &=& h\, F_2(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ k_4 &=& h\, F_1(x_i + k_3, v_i + m_3, t_i + h) \\ m_4 &=& h\, F_2(x_i + k_3, v_i + m_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4)\\ v_{i+1} &=& v_i + \frac{1}{6} (m_1 + 2 m_2 + 2 m_3 + m_4) \end{eqnarray}簡単な減衰＋強制振動の方程式
$$ \frac{d^2 x}{dt^2} = -x - a \frac{dx}{dt} + b \cos(t) $$を，
して解きます。
もう初めから組込関数 rk を使います。計算の精度は rk の引数である刻み幅 h に依存しますから，h の大きさを変えて調べてみます。
kill(x, v, t, a, b, N, n)$
F1(x, v, t):= v$
F2(x, v, t, a, b):= -x - a*v + b*cos(t)$
t0: 0 $
x0: 3 $
v0: 0$
t1: 20 $
N: 200$
h: float((t1 - t0) * 1.0/N)$
data00b00: rk([F1(x,v,t),F2(x,v,t,0,0)], [x,v], [x0,v0], [t, t0, t1, h])$
printf(true, "~10f ~10f~%", h, data00b00[length(data00b00)][2])$
h:h/10$
data00b00: rk([F1(x,v,t),F2(x,v,t,0,0)], [x,v], [x0,v0], [t, t0, t1, h])$
printf(true, "~10f ~10f~%", h, data00b00[length(data00b00)][2])$
h:h/10$
data00b00: rk([F1(x,v,t),F2(x,v,t,0,0)], [x,v], [x0,v0], [t, t0, t1, h])$
printf(true, "~10f ~10f~%", h, data00b00[length(data00b00)][2])$
以上の結果から，今回は h = 0.01 で小数点以下8桁くらいまでの精度が出ていると推測されます。以下ではこの h で計算を続けます。
h: 0.01$
data00b00: rk([F1(x,v,t),F2(x,v,t,0,0)], [x,v], [x0,v0], [t, t0, t1, h])$
data05b02: rk([F1(x,v,t),F2(x,v,t,0.5,0.2)], [x,v], [x0,v0], [t, t0, t1, h])$
組込関数 rk が返す計算結果のリストから，横軸・縦軸に表示するデータを makelist で取り出し，2つのデータのグラフを重ねて表示します。
plot2d([[discrete, makelist([c[1],c[2]], c, data00b00)], /* 3コラムあるデータのうち，1コラム目と2コラム目を */
[discrete, makelist([c[1],c[2]], c, data05b02)]], /* 使って plot2d する例 */
grid2d,
[legend, "a = 0 , b = 0", "a = 0.5, b = 0.2"])$
実は Maxima では，1階および2階常微分方程式を解析的に解く関数 ode2() が組み込まれています。
以下では， $$ \frac{d^2 x}{dt^2} = -x - a \frac{dx}{dt} + b \cos(t) $$ を，$a = 0.5, b = 0.2$ の場合に解析的に解く例を示します。
解くべき微分方程式の記述例です。微分 diff() の前に ' をつけて書きます。
a:1/2$
b:1/5$
eq: 'diff(x, t, 2) = - x - a*'diff(x,t) + b*cos(t);
ode2() 関数を使って常微分方程式を解く書式です。
ode2(eq, x, t);
上記の解は，2つの積分定数 $\%k_1, \%k_2$ を含む一般解です。
2階常微分方程式の初期条件は ic2() 関数で与えます。以下では，初期条件として $t = 0$ で $x = 3, \frac{dx}{dt} = 0$ を設定しています。
ic2(%, t = 0, x = 3, 'diff(x, t) = 0);
この解を関数 x(t) として定義し，plot2d() してみます。
define(x(t), rhs(%))$
plot2d(x(t), [t, 0, 20], grid2d)$
ルンゲ・クッタ法 rk() で求めた数値解と一緒に図示します。数値解は2001個と多いので，間引いたリスト rkdat を作り，描画スタイルを points にしてグラフにします。
length(data05b02);
rkdat: makelist(data05b02[i], i, 1, length(data05b02), 100)$
length(rkdat);
plot2d([x(t),[discrete, makelist([c[1],c[2]], c, rkdat)]], [t, 0, 20],
grid2d,
[style, lines, [points, 1]],
[legend, "解析解", "数値解"])$
念のために，$t=20$ での数値解と解析解を比べてみます。
printf(true, "~16h~%",data05b02[2001][2])$
printf(true, "~16h~%", x(20))$