Return to Maxima でコンピュータ演習

Maxima で微分積分・方程式の解

Maxima の基本から，関数，微分・積分，方程式の解まで。

JupyterHub での計算の実行と行末記号

  1. Maxima ではセルに実行文を書き，
  2. 行末にセミコロン ; を入力後，
  3. Shift キーを押しながら Enter キー（Return キー）を押すと実行します。
  4. SHift + Enter のかわりに，メニューの「▶︎ Run」をクリックしてもよいです。
  5. 行末にセミコロン ; の代わりに $ を入力後，実行すると，結果が表示されません。

まずは簡単なかけ算 123×456 を計算してみます。

In [1]:
123 * 456;
Out[1]:
\[\tag{${\it \%o}_{1}$}56088\]

2 の 100 乗 ($2^{100}$) のような大きい数でも，厳密な結果を出力します。（** は累乗を表します。^ でもいいですが，他の言語では ** を使うことが多いので，こちらに統一します。）

In [2]:
2**100;
Out[2]:
\[\tag{${\it \%o}_{2}$}1267650600228229401496703205376\]

float 関数を使って浮動小数点表示で近似値を出力させることもできます。

In [3]:
float(2**100);
Out[3]:
\[\tag{${\it \%o}_{3}$}1.26765060022823 \times 10^{+30}\]

表記などの約束事

Maxima の実行に際して，いくつか表記についての約束事があります。代表的なものを以下に示します。

コメント(実行に影響しない注釈)を 入れたいときは，/* ... */ でコメント部分を囲みます。

In [4]:
1 + 2 + 3;
Out[4]:
\[\tag{${\it \%o}_{4}$}6\]
In [5]:
1 + 2 /* + 3 */;
Out[5]:
\[\tag{${\it \%o}_{5}$}3\]

加法 +，減法 -，乗法 *，除法 / における優先順位は数学と同じです。優先して計算したい箇所は ( ) で 囲みます。

In [6]:
/* 暗算でもできる！*/
1 + 2 * 3;
Out[6]:
\[\tag{${\it \%o}_{6}$}7\]
In [7]:
(1 + 2) * 3;
Out[7]:
\[\tag{${\it \%o}_{7}$}9\]

ただし，大カッコ [ ] や中カッコ { } は優先順位の指定には用いません。

代入

変数に式や値を代入するときは : (コロン)を使います。

変数 a に，$100 + 3/21$ を代入して，
a を使って計算を続けます.

In [8]:
a: 100 + 3/21;
Out[8]:
\[\tag{${\it \%o}_{8}$}\frac{701}{7}\]
In [9]:
(a - 100) * 7;
Out[9]:
\[\tag{${\it \%o}_{9}$}1\]
In [10]:
a;
Out[10]:
\[\tag{${\it \%o}_{10}$}\frac{701}{7}\]

関数の定義

Maxima では，関数の定義には := を使います。$f(x) := x^2$

In [11]:
f(x):= x**2;
Out[11]:
\[\tag{${\it \%o}_{11}$}f\left(x\right):=x^2\]

$f(4) = 4^2 = 16$

In [12]:
f(4);
Out[12]:
\[\tag{${\it \%o}_{12}$}16\]

実行結果の非表示

文末に $ をつけると実行結果の出力を省略します。

In [13]:
eq: x^2 + 3*x + 2$

変数 eq を出力させると…
確かに式 $x^2 +3x+2$ ですね.

In [14]:
eq;
Out[14]:
\[\tag{${\it \%o}_{14}$}x^2+3\,x+2\]

前の出力結果の参照

直前の出力結果は % と表します。これを参照してさらに計算することができます。また，以前の結果を参照する場合は，例えば %o14 のように行番号を指定することもできます。

expand 関数を使って展開します。

In [15]:
expand((x - 2)**4);
Out[15]:
\[\tag{${\it \%o}_{15}$}x^4-8\,x^3+24\,x^2-32\,x+16\]

直前の結果を factor 関数で因数分解します。

In [16]:
factor(%);
Out[16]:
\[\tag{${\it \%o}_{16}$}\left(x-2\right)^4\]
In [17]:
/* 行番号を指定して以前の結果を表示させてみてください。*/

factor(%o15);
Out[17]:
\[\tag{${\it \%o}_{17}$}\left(x-2\right)^4\]

数の計算・基本的な関数

基本的な四則演算 (+ - * /) 以外によく使われる演算記号を示します。

4 の階乗の計算をします. 4! = 4*3*2*1 = 24 です。

In [18]:
4!;
Out[18]:
\[\tag{${\it \%o}_{18}$}24\]

!! は，一つおきの階乗の計算をします。5!! = 5*3*1 です。

In [19]:
5!!;
Out[19]:
\[\tag{${\it \%o}_{19}$}15\]
In [20]:
6!!;
Out[20]:
\[\tag{${\it \%o}_{20}$}48\]
In [21]:
6 * 4 * 2 ;
Out[21]:
\[\tag{${\it \%o}_{21}$}48\]

練習 1-1

(練習) 3, 2, 1 でつくる最大の数

数字の 3 と 2 と 1 を 1 個ずつ使った演算でつくることができる最大の整数は?

ヒント:

(1 + 2) × 3 = 9 じゃ全然だめ。そのまま並べると 321 ですが，累乗を使うともっと大きい数をつくることができます。
$21^3$ とか，$31^2$ とか，$2^{31}$ とか…

In [ ]:
基本的な定数
Maxima の基本的な定数を示します。
　　円周率 $\pi$　　　　無限大 $\infty$　　　　自然対数の底 $e$　　　　虚数単位 $i$　　
%piinf%e%i
無理数 $e$，$\pi$ と虚数単位 $i$ と整数 $−1$ の間のすてきな関係，オイラーの等式 $$e^{i\pi} = -1$$
In [22]:
%e**(%i * %pi);

Out[22]:
\[\tag{${\it \%o}_{22}$}-1\]
参考：display() を使って表示させる例。評価する前の式を左辺に，評価後の値を右辺にした式を表示します。
In [23]:
display(%e**(%i * %pi))$

\[e^{i\,\pi}=-1\]
浮動小数点で近似値を表示するには float を使います。
In [24]:
%pi; 
float(%pi);

Out[24]:
\[\tag{${\it \%o}_{24}$}\pi\]
Out[24]:
\[\tag{${\it \%o}_{25}$}3.141592653589793\]
基本的な関数
基本的な関数は以下のように表記されます。
平方根　指数関数　自然対数　絶対値　三角関数　逆三角関数　双曲線関数
sqrtexplogabssin
 cos
 tan		asin
 acos
 atan		sinh
 cosh
 tanh
うっかり $\sqrt{x^2} = x$ としないように。 $\sqrt{x^2} = |x|$ です。
In [25]:
sqrt(x**2);

Out[25]:
\[\tag{${\it \%o}_{26}$}\left| x\right| \]
ピタゴラスの定理をみたす整数の例。
 $\sqrt{3^2 + 4^2} = 5$ ですね。
In [26]:
sqrt(3**2 + 4**2);

Out[26]:
\[\tag{${\it \%o}_{27}$}5\]
$\exp(1) = e$ のこと。
$e^{\log x}$ や $\log(e^2)$ の計算。
In [27]:
display(exp(log(x)))$
display(log(%e**2))$

\[\exp \log x=x\]
\[\log e^2=2\]
三角関数の値・引数の単位
三角関数の値を数値で求めたいときには，以下のようにします。
$\displaystyle \sin\frac{\pi}{3}$ の値は…
In [28]:
sin(%pi/3);

Out[28]:
\[\tag{${\it \%o}_{30}$}\frac{\sqrt{3}}{2}\]
float 関数を使って，直前の結果を
 浮動小数点表示します。
In [29]:
float(%);

Out[29]:
\[\tag{${\it \%o}_{31}$}0.8660254037844386\]
変数 pi に浮動小数点表示の $\pi$ の値を代入します。 浮動小数点表示で数値を入れると， sin も浮動小数点表示で答えを返します。
In [30]:
pi: float(%pi);
sin(pi/3);

Out[30]:
\[\tag{${\it \%o}_{32}$}3.141592653589793\]
Out[30]:
\[\tag{${\it \%o}_{33}$}0.8660254037844386\]
三角関数の引数はラジアンです。度を使って求めたいときは以下のようにします。
度をラジアンに変換する関数 degrad を定義します。
degrad(60) を引数にして sin の値を求めます。
In [31]:
degrad(x):= x/180*%pi$
sin(degrad(60));

Out[31]:
\[\tag{${\it \%o}_{35}$}\frac{\sqrt{3}}{2}\]
In [32]:
sin(degrad(30));

Out[32]:
\[\tag{${\it \%o}_{36}$}\frac{1}{2}\]
練習 1-2
（練習）屋根勾配
家の屋根の勾配は角度ではなく，伝統的に 3 寸勾配，4 寸勾配のように水平方向に 1 尺(10 寸)に対して垂直方向に何寸立ち上がるかで示します。
図のログハウスの屋根は 8 寸勾配です。
 傾斜角度にすると何度ですか?
以下に各寸勾配の角度の表があります。10 寸勾配まで表
 を完成させましょう。
勾配　　1寸勾配　　　2寸勾配　　　3寸勾配　　　…
atan(0.1)atan(0.2)atan(0.3)…
rad0.09970.19740.2915…
度5.7106°11.3099°16.6992°…
In [ ]:
練習 1-3
(練習)加法定理・三倍角の公式
  1. 三角関数及び双曲線関数の加法定理を確認しなさい。
  2. $\sin 3x$ を $\sin x$ だけを使って表しなさい。
  3. $\tan 3x$ を $\tan x$ だけを使って表しなさい.
In [33]:
/* 1. のヒント */
/* trigexpand で展開します。*/
sin(x + y) = trigexpand(sin(x + y))$

In [34]:
/* 2. や 3. のヒント */
sin(3*x) = trigexpand(sin(3*x))$
trigsimp(%)$

$\cos(x)^2$ を $1-\sin(x)^2$ で置き換える例です。
In [35]:
ratsubst(1-sin(x)^2, cos(x)^2, %)$

In [ ]:
微分・積分・数値積分
微分・積分を行う関数の表記を以下に示します。
操作 　表記例
微分diff(f(x), x)
n 階微分diff(f(x), x, n)
不定積分integrate(f(x), x)
定積分integrate(f(x), x, a, b)
数値積分romberg(f(x), x, a, b), quad_qags(f(x), x, a, b)
微分・積分
理工系の数学Bで Maxima を使った微分・積分についてやっています。以下を参照のこと。
数値積分：romberg()
では，解析的に積分できない場合はどうするか？
実際，大学で習う「初等関数」およびその合成関数は全て微分できますが，では初等関数およびその合成関数（初等関数を使って書かれる関数）は全て積分できるかというと，そうではない。
大学の授業，たとえば理工系の数学Bなどでは，解析的に積分できるような例題だけを扱いますが，実際には，以下の例のように $\sin(\sin x)$ などのような単純な例でも，解析的には積分できません。
解析的に積分できない場合は，数値積分によって近似値を求めます。Maxima には数値積分を行う関数があります。
$\displaystyle \int_0^{\pi} \sin(\sin x)\ dx$ の計算をします。
In [36]:
integrate(sin(sin(x)), x, 0, %pi);

Out[36]:
\[\tag{${\it \%o}_{41}$}\int_{0}^{\pi}{\sin \sin x\;dx}\]
解析的に解けないと，そのままで返します。
romberg 関数を使うと，数値的に積分した近似解を表示します。
In [37]:
r1: romberg(sin(sin(x)), x, 0, %pi);

Out[37]:
\[\tag{${\it \%o}_{42}$}1.786487487137068\]
数値積分 romberg の精度は，大域変数 rombergtol で決まります。
rombergtol の値を小さくすることで精度をあげることができます。
デフォルトの rombergtol の値は…
In [38]:
rombergtol;

Out[38]:
\[\tag{${\it \%o}_{43}$}1.0 \times 10^{-4}\]
rombergtol を $10^{-10}$ にしてもう一度同じ数値積分を行うと…
In [39]:
rombergtol: 1.e-10$
r2: romberg(sin(sin(x)), x, 0, %pi);

Out[39]:
\[\tag{${\it \%o}_{45}$}1.78648748195006\]
rombergtol を $10^{-12}$ にしてもう一度同じ数値積分を行うと…
In [40]:
rombergtol: 1.e-12$
r3: romberg(sin(sin(x)), x, 0, %pi);

Out[40]:
\[\tag{${\it \%o}_{47}$}1.786487481950052\]
以上のことから，数値積分は，あくまで近似的解を求めるものだが，どのくらいまで正しく，どのくらいの誤差があるかは，自分でコントロールできる。rombergtol はデフォルトでは $10^{-4}$ なので，もっと精度を求める場合は $10^{-10}$ とか $10^{-12}$ にしておけば，それなりの精度をもつ数値解が求められるということになる。
数値積分：quad_qags()
Maxima には数値積分を行う関数として quad_qags() もあります。romberg() に対する quad_qags() の利点については，以下のリンクに簡単にまとめています。
quad_qags() で $\displaystyle \int_0^{\pi} \sin(\sin x)\ dx$ の計算をしてみます。
In [41]:
quad_qags(sin(sin(x)), x, 0, %pi);

Out[41]:
\[\tag{${\it \%o}_{48}$}\left[ 1.786487481950052 , 2.468727406962818 \times 10^{-11} , 21 , 0 \right] \]
quad_qags() は4つの要素のリストを返します。
  • 積分の近似値
  • 見積もられた近似の絶対誤差
  • 被積分関数の評価数
  • エラーコード
デフォルトでは $10^{-11}$ 程度の誤差の数値積分値を返します。
数値積分の精度を設定し（デフォルトでは epsrel=1d-8 ），積分値のみを取り出す例です。
In [42]:
ans: 
  quad_qags(sin(sin(x)), x, 0, %pi, 
            'epsrel=1d-12)$
ans[1];

Out[42]:
\[\tag{${\it \%o}_{50}$}1.786487481950052\]
方程式
多項式方程式・連立方程式
Maxima は，式の計算だけでなく，代数方程式を解くことができます.例を以下に示します。
$x^2 +2x−2 = 0$ を $x$ について解き
 ます.
In [43]:
solve(x**2 + 2*x - 2 = 0, x);

Out[43]:
\[\tag{${\it \%o}_{51}$}\left[ x=-\sqrt{3}-1 , x=\sqrt{3}-1 \right] \]
連立方程式 $x^2 + y^2 = 5, \ x + y = 1$ を $x, y$ について解きます。
In [44]:
solve(
  [x**2 + y**2 = 5,  x + y = 1], 
  [x, y]
);

Out[44]:
\[\tag{${\it \%o}_{52}$}\left[ \left[ x=-1 , y=2 \right] , \left[ x=2 , y=-1 \right] \right] \]
3 次方程式 $x^3 = 8$ を解きます。
In [45]:
solve(x**3 = 8, x);

Out[45]:
\[\tag{${\it \%o}_{53}$}\left[ x=\sqrt{3}\,i-1 , x=-\sqrt{3}\,i-1 , x=2 \right] \]
確認のため，上記の1番目の答えを3乗して簡単化すると…
In [46]:
/* 直前の結果は [] で囲まれたリストになっています。*/
/* リストの1番目をとるには %[1] のようにします。*/
rhs(%[1])**3;

Out[46]:
\[\tag{${\it \%o}_{54}$}\left(\sqrt{3}\,i-1\right)^3\]
In [47]:
/* 展開します。*/
expand(%);

Out[47]:
\[\tag{${\it \%o}_{55}$}8\]
練習 1-4
（練習）鶴と亀と蟻
鶴と亀と蟻，個体数の合計は 10。足の数は全部で 34 本。蟻は亀より 1 匹少ない。鶴，亀，蟻，それぞれ何羽・何匹?
鶴と亀と蟻の個体数をそれぞれ $x, y, z$ として連立方程式をたてて solve します。蟻の足は何本かわかりますよね?
かつて線形代数の授業でこの問題を出したとき，真っ先に出た質問が「先生，蟻の足は何本ですか？」でした… orz
In [ ]:
方程式の数値解
多項式方程式でない場合には find_root 関数を使って，範囲を指定して解を探します。
関数 $g(x):= (x − 1)\,e^x$ の定義。
$g(x):= (x − 1)\,e^x = 1$ の解は solve などでは求めることができません。
In [48]:
g(x):= (x - 1)*exp(x);
solve(g(x) = 1, x);

Out[48]:
\[\tag{${\it \%o}_{56}$}g\left(x\right):=\left(x-1\right)\,\exp x\]
Out[48]:
\[\tag{${\it \%o}_{57}$}\left[ x=e^ {- x }\,\left(e^{x}+1\right) \right] \]
find_root を使って $0 \leq x \leq 2$ の 範囲で解を探してみます。 （なぜこの範囲に解があるとわかるのか，は以下で。）
出た答えを代入して確かめると…
In [49]:
find_root(g(x) = 1, x, 0, 2);
g(%);

Out[49]:
\[\tag{${\it \%o}_{58}$}1.278464542761074\]
Out[49]:
\[\tag{${\it \%o}_{59}$}1.0\]
$g(x)$ がどのへんで $1$ になるかは，以下のようにグラフを描いてみればわかります。Maxima によるグラフ作成については，別途解説します。
In [50]:
plot2d(g(x), [x, 0, 5], [y, -1, 1], [ylabel, "g(x)"])$

練習 1-5
（練習）関数の極大値
次の関数が極大値をとる $x$ の値およびその極大値を求めよ。
$$ 1. \ \ f(x) = \frac{x^3}{e^x – 1}, \qquad 2. \ \ g(x) = \frac{x^5}{e^x -1}$$
In [52]:
/* 1. のヒント */
/* f(x) を定義して... */
f(x):= x**3/(exp(x) - 1)$

/* f(x) を微分して... */
diff(f(x), x)$

/* 微分した直前の結果を関数 df(x) と定義。*/
define(df(x), %)$

In [53]:
/* 導関数 df(x) のグラフを描き，あたりをつける。*/
plot2d(df(x), [x, 0.1, 10])$

In [55]:
/* 導関数 df(x) がゼロになる x を 0 < x < 4 の間で探す＞*/
find_root(df(x) = 0, x, 2, 4)$

In [56]:
/* 極値をとる x の値が求まったので，極値は... */
f(%)$

In [57]:
/* 2. も同様にやっておく。*/

参考：プランクの法則
熱力学，気象学，天文学，宇宙物理学等で現れる「プランクの法則」。
温度 $T$ の黒体から放射される電磁波の放射強度は周波数を $\nu$ として以下のように与えられる。
$$I(\nu, T) = \frac{2h\nu^3}{c^2} \frac{1}{e^{h\nu/kT} – 1}$$
ここで，$\displaystyle x \equiv \frac{h\nu}{kT}$ とおくと
$$I(\nu, T) = \frac{2 (k T)^3}{h^2 c^2} \frac{x^3}{e^{x} – 1} \equiv C \frac{x^3}{e^{x} – 1}$$
となり，$f(x)$ の極大値を調べることは，プランクの法則で与えられる黒体からの放射強度の極大値を調べることにつながるんですよ。