Return to 参考:Maxima 編(旧版アーカイブ)

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)$ です。)

In [1]:
f(x):= sin(x)$

/* f(x) の数値微分 */
/* Maxima では数値微分する関数はないので自作する */
diff_f(x, h):= (f(x+h) - f(x-h))/(2.0*h)$

printf(true, "数値微分  ~f~%", diff_f(3, 1e-6))$
printf(true, "解析解   ~f~%", cos(3.0))$
数値微分  -0.9899924967304852
解析解   -0.9899924966004454

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

In [2]:
plot2d(cos(x) - diff_f(x, 1.0E-6), [x, -2*%pi, 2*%pi], [ylabel,""])$

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

In [3]:
plot2d(cos(x) - diff_f(x, 1.0E-7), [x, -2*%pi, 2*%pi], [ylabel,""])$
system("cp ~/.maxplot.svg ./maxsk02.svg")$

数値積分

romberg()

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

In [4]:
romberg(sin(sin(x)), x, 0, %pi);
Out[4]:
\[\tag{${\it \%o}_{9}$}1.786487487137068\]

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

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

デフォルトの rombergtol の値は…

In [5]:
rombergtol;
Out[5]:
\[\tag{${\it \%o}_{10}$}1.0 \times 10^{-4}\]

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

rombergtol: 1e-11$

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

In [6]:
for h in [1e-8, 1e-9, 1e-10, 1e-11, 1e-12, 1e-13] do(
  rombergtol: h,
  printf(true, "~17,15f~%", romberg(sin(sin(x)), x, 0, %pi))
)$
1.786487481932776
1.786487481950060
1.786487481950060
1.786487481950052
1.786487481950052
1.786487481950052

quad_qags()

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

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

In [7]:
quad_qags(sin(sin(x)), x, 0, %pi);
Out[7]:
\[\tag{${\it \%o}_{12}$}\left[ 1.786487481950052 , 2.468727406962818 \times 10^{-11} , 21 , 0 \right] \]

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

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

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

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

In [8]:
ans: 
  quad_qags(sin(sin(x)), x, 0, %pi, 'epsrel=1d-12)$
ans[1];
Out[8]:
\[\tag{${\it \%o}_{14}$}1.786487481950052\]

方程式の数値解 find_root()

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

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

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

In [9]:
f(x):= bessel_j(0, x)$
plot2d(f(x), [x, 0, 10], [legend, "J_0 ベッセル関数"], grid2d)$

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

In [10]:
find_root(f(x) = 0, x, 2, 3);
Out[10]:
\[\tag{${\it \%o}_{18}$}2.404825557695773\]
In [11]:
/* 本当にゼロになっているか,直前の答え % を入れてみます */
f(%);
Out[11]:
\[\tag{${\it \%o}_{19}$}-5.551115123125783 \times 10^{-17}\]

上の例で,% は直前の出力(計算結果)を表します。

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() があらかじめ組み込まれています。使い方は以下のとおりです。

In [12]:
/* dx/dt の右辺 */
f(x, t):= (1-x)*x$

/* 初期条件と計算範囲 */
t0: 0$
x0: 0.1$
t1: 10$

/* きざみ幅 h */
N: 100$
h: float((t1 - t0) * 1.0/N)$

results: rk(f(x, t), x, x0, [t, t0, t1, h])$
/* h と計算結果を表示 */
printf(true, "h = ~5,4f   x(t1) = ~f~%", h, results[length(results)][2])$
h = .1000   x(t1) = 0.9995915652218679

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

In [13]:
for N in [400, 1000, 4000, 10000] do(
  h: float((t1 - t0) * 1.0/N),
  results: rk(f(x, t), x, x0, [t, t0, t1, h]),
  printf(true, "h = ~7,5f   x(t1) = ~f~%", h, results[length(results)][2])
)$
h = 0.02500   x(t1) = 0.9995915675088635
h = 0.01000   x(t1) = 0.9995915675171755
h = 0.00250   x(t1) = 0.999591567517391
h = 0.00100   x(t1) = 0.9995915675173903

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

数値解をグラフに

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

In [14]:
plot2d([discrete, results], [y, 0, 1], grid2d)$

In [15]:
/* plot2d() のかわりに draw2d() を使う例 */
draw2d(
  grid = true,
  points_joined = true, point_type = dot,
  points(results)
)$

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 の大きさを変えて調べてみます。

In [16]:
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)$ 

for i:1 thru 4 do(
  data00b00: rk([F1(x, v, t), F2(x, v, t, 0, 0)], 
                [x, v], [x0, v0], [t, t0, t1, h]),
  printf(true, "h = ~6,4f   x(t1) = ~15f~%", 
                h, data00b00[length(data00b00)][2]),
  h: h/10
)$
h = 0.1000   x(t1) = 1.2242899713355
h = 0.0100   x(t1) = 1.2242461899877
h = 0.0010   x(t1) = 1.2242461854406
h = 0.0001   x(t1) = 1.2242461854401

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

In [17]:
h: 0.01$
/* a = 0, b = 0 */
data00b00: rk([F1(x, v, t), F2(x, v, t, 0, 0)], 
              [x, v], [x0, v0], [t, t0, t1, h])$
/* a = 0.5, b = 0.2 */
data05b02: rk([F1(x, v, t), F2(x, v, t, 0.5, 0.2)], 
              [x, v], [x0, v0], [t, t0, t1, h])$

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

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

In [18]:
/* 3コラムあるデータのうち,1コラム目と2コラム目を */
/* 使って plot2d する例 */
plot2d([[discrete, makelist([c[1],c[2]], c, data00b00)],  
        [discrete, makelist([c[1],c[2]], c, data05b02)]], 
       grid2d, [y, -3, 3.9], /* 凡例が曲線とカブるのを避ける */
       [legend, "a = 0 , b = 0", "a = 0.5, b = 0.2"])$