Maxima で数値解析

あちらの「Maxima で数値解析」では，Maxima の数値解析的な関数を利用するほかに，解析的な方法や自前で数値解析的コードを書く方法なども併せて紹介しているので，ここでは数値解析な関数の利用について抜粋して紹介する。

数値微分

解析的な微分の定義は（前方差分の極限として）
$$\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 を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。

\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}

例として，$f(x) = \sin(x)$ の $x=3$ における数値微分を求めてみます。（答えは $\cos(3)$ です。）

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

以下の例からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけでもなさそうです。

数値積分

romberg()

$\displaystyle \int_0^{\pi} \sin (\sin x) \, dx$ のように解析的に積分できない場合は，romberg() という組込関数で数値積分できます。

数値積分 romberg の精度は，大域変数 rombergtol で決まります。

rombergtol の値を小さくすることで精度をあげることができます。

デフォルトの rombergtol の値は…

rombergtol の値を変えて計算し，精度を確認します。以下の結果から，

rombergtol: 1e-11$

で小数点以下15桁程度の精度があることがわかります。

quad_qags()

Maxima には数値積分を行う関数として quad_qags() もあります。romberg() に対する quad_qags() の利点については，以下のリンクに簡単にまとめています。

• 数値積分 quad_qags()

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

quad_qags() は4つの要素のリストを返します。

• 積分の近似値
• 見積もられた近似の絶対誤差
• 被積分関数の評価数
• エラーコード

デフォルトでは $10^{-11}$ 程度の誤差の数値積分値を返します。

数値積分の精度を設定し，積分値のみを取り出す例です。

方程式の数値解 find_root()

$f(x) = 0$ の解を数値的に求めるには， find_root() を使います。

まず，解を求める関数 $f(x)$ を定義し，$f(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。

以下では例として組み込み関数の1つであるベッセル関数 $J_0(x)$ bessel_j(0,x) のゼロ点を求めています。

以下は，$2 < x < 3$ の範囲で $f(x) = 0$ の解を数値的に求める例です。

上の例で，% は直前の出力（計算結果）を表します。

1階常微分方程式の数値解 rk()

簡単なロジスティック方程式
$$ \frac{dx}{dt} = f(x, t) = (1-x)\,x$$
を，初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 0.1$，$t = t_0 = 0$ から $t_1 = 10$ までを数値的に解きます。

Maxima にはルンゲ・クッタ法によって常微分方程式を解く関数 rk() があらかじめ組み込まれています。使い方は以下のとおりです。

Maxima の rk() では h の大きさが直接精度に関係します。N の値を変えて h の値を変え，計算精度を確認します。

以上の結果から，例えば小数点以下14桁程度の精度がどうしても欲しいということであれば，h = 0.001 程度にする必要があることがわかります。

数値解をグラフに

数値計算の結果をグラフにしてみます。

2階常微分方程式の数値解 rk()

次のような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}

簡単な減衰＋強制振動の方程式

$$ \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 $ 分割

して解きます。

計算の精度は rk の引数である刻み幅 h に依存しますから，h の大きさを変えて調べてみます。

以上の結果から，今回は h = 0.01 で小数点以下8桁くらいまでの精度が出ていると推測されます。以下ではこの h で計算を続けます。

2つの数値解を重ねてグラフに

組込関数 rk が返す計算結果のリストから，横軸・縦軸に表示するデータを makelist で取り出し，2つのデータのグラフを重ねて表示します。