Maxima-Jupyter による数式処理とグラフ作成

この Notebook では,「wxMaxima による数式処理とグラフ作成」のテキストの内容を,Jupyter Notebook 上で Maxima-Jupyter を使って説明しています。

セクションの構成は,wxMaxima 版のテキストに準じています。


葛西 真寿(KASAI Masumi)           
弘前大学 大学院理工学研究科        
情報連携統括本部 情報基盤センター

Maxima-Jupyter による数式処理

基本操作

Maxima-Jupyter の起動

ターミナル等で

jupyter notebook (または jupyter-notebook)

右上の「New」から 「Maxima」を選ぶ。

計算の実行と改行

計算の実行は,数値や式などを入力して最後にセミコロン ; を入力後,Shift キーを押しならが Enter キー(または Return キー)を押します。

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

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

2 の 100 乗 (2100) のような大きい数でも,厳密な結果を出力します。(^ は累乗を表します。)

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

Enter キーのみでは改行されるだけで計算が実行されません。Shift+Enter キーが入力されるまで,Maxima は式の入力が継続していると判断します。

例として,以下のような積分を計算してみます。Shift キーを忘れて Enter キーのみを押すと単に改行されます。慌てずにセミコロン ; を入力後, Shift キーを押しながら Enter キーを押すと実行します。

答えは $\int 3x^2\ dx = x^3$ です.

In [4]:
integrate(3*x^2, x)


;
Out[4]:
\[\tag{${\it \%o}_{4}$}x^3\]

計算結果の保存と読み込み

Jupyter Notebook (の Maxima-Jupyter)で計算した結果を保存するには,「File」メニューから「Save as ... 」でファイル名をつけて保存します。

「File」メニューの真下に,なにやら四角で右上隅がかけている形のアイコンがありますが(これが「フロッピーディスク」というものであることは,20世紀の昔を知る人にしかわからないと思います),このアイコンをクリックしても保存されます。

また,「File」メニューから「Open... 」を選び,計算結果が保存された Notebook をまた利用することもできます。

表記などの約束事

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

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

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

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

In [7]:
1 + 2 * 3;
Out[7]:
\[\tag{${\it \%o}_{7}$}7\]
In [8]:
(1 + 2) * 3;
Out[8]:
\[\tag{${\it \%o}_{8}$}9\]

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

代入

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

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

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

変数 eq に式 $x^2 + 3 x + 2$ を代入します。

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

eqfactor 関数で因数分解します。

In [12]:
factor(eq);
Out[12]:
\[\tag{${\it \%o}_{12}$}\left(x+1\right)\,\left(x+2\right)\]

関数の定義

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

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

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

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

$f(3\sqrt{A}) = (3\sqrt{A})^2 = 9 A$ です。

In [15]:
f(sqrt(A)*3);
Out[15]:
\[\tag{${\it \%o}_{15}$}9\,A\]

$\frac{d}{dx}(x \cos x)$ を計算します。

In [16]:
diff(x*cos(x), x);
Out[16]:
\[\tag{${\it \%o}_{16}$}\cos x-x\,\sin x\]

以下の例のように f(x):= % のような 仕方で(直前の)計算結果を参照して関数を定義することはできません。

In [17]:
g(x):= %;
Out[17]:
\[\tag{${\it \%o}_{17}$}g\left(x\right):={\it \%}\]
In [18]:
g(x);
Out[18]:
\[\tag{${\it \%o}_{18}$}g\left(x\right):={\it \%}\]

このような場合には define を使います。%o16 は 16 番目の出力結果(今の場合は $\cos(x) − x \sin(x)$ です。)

In [19]:
define(g(x), %o16);
Out[19]:
\[\tag{${\it \%o}_{19}$}g\left(x\right):=\cos x-x\,\sin x\]
In [20]:
integrate(g(x), x);
Out[20]:
\[\tag{${\it \%o}_{20}$}x\,\cos x\]

実行結果の非表示

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

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

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

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

2 次方程式 $x^2 +3x+2 = 0$ を solve を使って解いていますが,結果は出力されません。$ ではなく ; を文末 につけると結果が出力されます。

In [23]:
solve(eq = 0, x)$
In [24]:
solve(eq = 0, x);
Out[24]:
\[\tag{${\it \%o}_{24}$}\left[ x=-2 , x=-1 \right] \]

前の出力結果の参照

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

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

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

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

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

行番号 %o26 の結果を展開します。

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

一時的代入

ev 関数を使うことで,一時的に式などに値を代入する事ができます。

変数 c2*x を代入します。

In [28]:
c: 2*x;
Out[28]:
\[\tag{${\it \%o}_{28}$}2\,x\]

x=3 を一時的に代入して c の値を評 価します。

In [29]:
ev(c, x=3);
Out[29]:
\[\tag{${\it \%o}_{29}$}6\]

c の定義そのものは変わりません。

In [30]:
c;
Out[30]:
\[\tag{${\it \%o}_{30}$}2\,x\]

リスト成分や右辺の取り出し

方程式の解などは,リスト形式 ([ ] で囲まれた形)で出力します。 リスト成分は [ ] で取 り出すことができます。

関数 f(x) を定義します。

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

solve で方程式 f(x) = 0 を解きます。

In [32]:
sol: solve( f(x) = 0, x );
Out[32]:
\[\tag{${\it \%o}_{32}$}\left[ x=-2 , x=-1 \right] \]

リストとして表されている 2 個の解のうち,1 番目を sol[1] として, 2 番目を sol[2] としてそれぞれ取り出します。

In [33]:
sol[1];
sol[2];
Out[33]:
\[\tag{${\it \%o}_{33}$}x=-2\]
Out[33]:
\[\tag{${\it \%o}_{34}$}x=-1\]

rhs は右辺 (right hand side) を取り出す関数です。また,左辺を取り出す関数は lhs (left hand side) です。

x = -2 の右辺を取り出して ans1 に代入します。

二つめの解も同様に ans2 に代入します。

In [34]:
ans1: rhs(sol[1]);
ans2: rhs(sol[2]);
Out[34]:
\[\tag{${\it \%o}_{35}$}-2\]
Out[34]:
\[\tag{${\it \%o}_{36}$}-1\]

関数 f(x) に解を入れると確かに 0 になります。

In [35]:
f(ans1);
Out[35]:
\[\tag{${\it \%o}_{37}$}0\]
In [36]:
f(ans2);
Out[36]:
\[\tag{${\it \%o}_{38}$}0\]

以下は見慣れない書式ですが, ev(f(x), sol[1]) の簡略形です。

In [37]:
f(x), sol[1];
Out[37]:
\[\tag{${\it \%o}_{39}$}0\]

数の計算・基本的な関数

基本的な演算

1 から 10 までの和を求めてみます。

In [38]:
1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10;
Out[38]:
\[\tag{${\it \%o}_{40}$}55\]

同様の計算は公式 $\sum_{i = 1}^{n} = n (n + 1)/2$ を 使ってもできます。

In [39]:
n: 10$
n*(n + 1)/2;
Out[39]:
\[\tag{${\it \%o}_{42}$}55\]
In [40]:
kill(n)$
ev(n*(n + 1)/2, n = 10);
Out[40]:
\[\tag{${\it \%o}_{44}$}55\]

以下は,ev(n*(n + 1)/2, n = 10) の簡略形でしたね。

In [41]:
n*(n+1)/2, n = 10;
Out[41]:
\[\tag{${\it \%o}_{45}$}55\]

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

$2^4$ の計算。累乗は ^ で表します。^ のかわりに ** でもよいようです。

In [42]:
2^4;
Out[42]:
\[\tag{${\it \%o}_{46}$}16\]
In [43]:
2**4;
Out[43]:
\[\tag{${\it \%o}_{47}$}16\]

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

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

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

In [45]:
5!!;
Out[45]:
\[\tag{${\it \%o}_{49}$}15\]

練習

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

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

ヒント:

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

基本的な定数

Maxima の基本的な定数を示します。

円周率 $\pi$ 無限大 $\infty$ 自然対数の底 $e$ 虚数単位 $i$
%pi inf %e %i

無理数 $e$,$\pi$ と虚数単位 $i$ と整数 $−1$ の間のすてきな関係,オイラーの等式 $$e^{i\pi} = -1$$

In [46]:
%e^(%i * %pi);
Out[46]:
\[\tag{${\it \%o}_{50}$}-1\]

浮動小数点で近似値を表示するには floatbfloat を使います。

In [47]:
%pi; 
float(%pi);
bfloat(%pi);
Out[47]:
\[\tag{${\it \%o}_{51}$}\pi\]
Out[47]:
\[\tag{${\it \%o}_{52}$}3.141592653589793\]
Out[47]:
\[\tag{${\it \%o}_{53}$}3.141592653589793_B \times 10^{0}\]

float は 16 桁の値を返します.それ以上の桁数を見たいときには,bfloat を使います。bfloat の表示桁数は fpprec で指定します。

有効数字 32 桁で円周率 $\pi$ を表示す る例です。

In [48]:
fpprec: 32$ bfloat(%pi);
Out[48]:
\[\tag{${\it \%o}_{55}$}3.1415926535897932384626433832795_B \times 10^{0}\]

基本的な関数

基本的な関数は以下のように表記されます。

平方根 指数関数 自然対数 絶対値 三角関数 逆関数 双曲線関数
sqrt exp log abs sin cos tan asin acos atan sinh cosh tanh

うっかり $\sqrt{x^2} = x$ としないように。 $\sqrt{x^2} = |x|$ です。

In [49]:
sqrt(x^2);
Out[49]:
\[\tag{${\it \%o}_{56}$}\left| x\right| \]

ピタゴラスの定理をみたす整数の例。 $\sqrt{3^2 + 4^2} = 5$ ですね。

In [50]:
sqrt(3^2 + 4^2);
Out[50]:
\[\tag{${\it \%o}_{57}$}5\]

$\exp(1) = e$ のこと。

$e^{\log x}$ や $\log(e^2)$ の計算。

In [51]:
exp(log(x));
log(%e^2);
Out[51]:
\[\tag{${\it \%o}_{58}$}x\]
Out[51]:
\[\tag{${\it \%o}_{59}$}2\]

練習

(練習)ピタゴラス数

$a^2 + b^2 = c^2$ を満たす正の整数の組をピタゴラス数といいます。ピタゴラス数は二つの任意の正の整数 $m, n$ から以下のように何組でもつくることができます。

$$a=|m^2 −n^2|, \ b=2mn, \ c= \sqrt{a^2 +b^2} $$

このとき,$c$ は必ず整数になることを示しなさい。

試しにやってみます。

In [52]:
kill(m, n)$ /* m, n を初期化 */
a: abs(m^2 - n^2);
b: 2 * m * n;
c: sqrt(a^2 + b^2);
Out[52]:
\[\tag{${\it \%o}_{61}$}\left| n^2-m^2\right| \]
Out[52]:
\[\tag{${\it \%o}_{62}$}2\,m\,n\]
Out[52]:
\[\tag{${\it \%o}_{63}$}\sqrt{\left(n^2-m^2\right)^2+4\,m^2\,n^2}\]

直接証明することも簡単ですが,疑似乱数を使って,ちょっとした数値実験をしてみます。

random(100)0 から 99 までの整 数疑似乱数をつくる関数です。

In [53]:
[a, b, c], m = random(100), n = random(100);
[a, b, c], m = random(100), n = random(100);
[a, b, c], m = random(100), n = random(100);
Out[53]:
\[\tag{${\it \%o}_{64}$}\left[ 140 , 48 , 148 \right] \]
Out[53]:
\[\tag{${\it \%o}_{65}$}\left[ 6069 , 5780 , 8381 \right] \]
Out[53]:
\[\tag{${\it \%o}_{66}$}\left[ 8265 , 728 , 8297 \right] \]

$c = \sqrt{a^2 + b^2}$ は確かに整数になっています。証明も簡単ですね。

練習

(練習)星の明るさと等級

夜空に見える星の明るさ $L$ は等級 $m$ で表され,以下のような関係があります。

$m = −2.5 \log_{10} L + \mbox{const.}\quad$ ($\mbox{const.}$ は定数のこと)

1 等星と 6 等星では,どちらが何倍明るいですか?

load(log10)$ とすると,常用対数 $\log_{10}$ が使えます.

$$\log_{10} 1000 = \log_{10} 10^3 = 3$$

ですね.

In [54]:
load(log10)$
log10(1000);
Out[54]:
\[\tag{${\it \%o}_{68}$}3\]

三角関数の値・引数の単位

三角関数の値を数値で求めたいときには,以下のようにします。

$\displaystyle \sin\frac{\pi}{3}$ の値は...

In [55]:
sin(%pi/3);
Out[55]:
\[\tag{${\it \%o}_{69}$}\frac{\sqrt{3}}{2}\]

float 関数を使って,直前の結果を 浮動小数点表示します。

In [56]:
float(%);
Out[56]:
\[\tag{${\it \%o}_{70}$}0.8660254037844386\]

変数 pi に浮動小数点表示の $\pi$ の値を代入します。 浮動小数点表示で数値を入れると, sin も浮動小数点表示で答えを返します。

In [57]:
pi: float(%pi);
sin(pi/3);
Out[57]:
\[\tag{${\it \%o}_{71}$}3.141592653589793\]
Out[57]:
\[\tag{${\it \%o}_{72}$}0.8660254037844386\]

三角関数の引数はラジアンです。度を使って求めたいときは以下のようにします。

度をラジアンに変換する関数 degrad を定義します。

degrad(60) を引数にして sin の値を求めます。

In [58]:
degrad(x):= x/180*%pi$
sin(degrad(60));
Out[58]:
\[\tag{${\it \%o}_{74}$}\frac{\sqrt{3}}{2}\]

練習

(練習)屋根勾配

家の屋根の勾配は角度ではなく,伝統的に 3 寸勾配,4 寸勾配のように水平方向に 1 尺(10 寸)に対して垂直方向に何寸立ち上がるかで示します。

図のログハウスの屋根は 8 寸勾配です。 傾斜角度にすると何度ですか?

以下に各寸勾配の角度の表があります。10 寸勾配まで表 を完成させましょう。

勾配 1寸勾配 2寸勾配 3寸勾配 ...
atan(0.1) atan(0.2) atan(0.3) ...
rad 0.0997 0.1974 0.2915 ...
5.7106° 11.3099° 16.6992° ...

数列の和

数列の和を計算する関数 sum の書式は,sum(一般項 a_i,i,はじめの i,最後の i) です。

$\displaystyle \sum_{i=1}^{100} i^2$ の計算をします。

In [59]:
sum(i^2, i, 1, 10);
Out[59]:
\[\tag{${\it \%o}_{75}$}385\]

一般の $n$ までの和は $\displaystyle \sum_{i=1}^{n} i^2$ ですが...

In [60]:
kill(n)$
sum(i^2, i, 1, n);
Out[60]:
\[\tag{${\it \%o}_{77}$}\sum_{i=1}^{n}{i^2}\]

sum 関数は答えを出してくれません。nusum は $n$ までの和を知っています。

In [61]:
nusum(i^2, i, 1, n);
Out[61]:
\[\tag{${\it \%o}_{78}$}\frac{n\,\left(n+1\right)\,\left(2\,n+1\right)}{6}\]

練習

(練習)数列の和

  1. 初項が $a_1$,公差 $d$ の等差数列の第 $1$ 項から第 $n$ 項までの和を求めよ。
  2. 初項が $b_1$,公比が $r$ の等比数列について,第 $n$ 項までの和を求めよ。
  3. 初項が $3$,公比が $2$ の等比数列について,第 10 項までの和を求めよ。
  4. 無限等比数列の和 $\displaystyle \sum_{i=1}^{\infty} b_1\,r^{i-1}$ を求めよ。ただし,$0 < r < 1$。
In [62]:
/* 1. のヒント */
a(a1, d, i):= a1 + d*(i-1);
nusum(a(a1, d, i), i, 1, n)$
Out[62]:
\[\tag{${\it \%o}_{79}$}a\left(a_{1} , d , i\right):=a_{1}+d\,\left(i-1\right)\]
In [63]:
/* 2. と 3. のヒント */
kill(b, b1, r, i)$

b(b1, r, i):= b1 * r^(i-1);
nusum(b(b1, r, i), i, 1, n)$
%, b1=3, r=2, n=10$
Out[63]:
\[\tag{${\it \%o}_{82}$}b\left(b_{1} , r , i\right):=b_{1}\,r^{i-1}\]
In [64]:
/* 4. のヒント */
/* r の範囲を宣言して... */
assume( 0 < r, r < 1);
nusum(b(b1, r, i), i, 1, n)$
limit(%, n, inf)$
Out[64]:
\[\tag{${\it \%o}_{85}$}\left[ r>0 , r<1 \right] \]

関数のテイラー展開

$\sin x$ を$x = 0$ のまわりで $x^5$ までテイラー展開する例。

In [65]:
taylor(sin(x), x, 0, 5);
Out[65]:
\[\tag{${\it \%o}_{88}$}\frac{x^5}{120}-\frac{x^3}{6}+x\]

$e^{ix}$($i$ は虚数単位)のテイラー展開。

In [66]:
ty: taylor(exp(%i*x), x, 0, 5);
Out[66]:
\[\tag{${\it \%o}_{89}$}\frac{i\,x^5}{120}+\frac{x^4}{24}-\frac{i\,x^3}{6}-\frac{x^2}{2}+i\,x+1\]

上記の虚部を imagpart,実部を realpart でとります。

In [67]:
imagpart(ty);
Out[67]:
\[\tag{${\it \%o}_{90}$}\frac{x^5}{120}-\frac{x^3}{6}+x\]
In [68]:
realpart(ty);
Out[68]:
\[\tag{${\it \%o}_{91}$}\frac{x^4}{24}-\frac{x^2}{2}+1\]

$\cos x, \tan x, \log (1 + x)$ 等のテイラー展開もやってみましょう。

In [69]:
taylor(log(1 + x), x, 0, 5);
Out[69]:
\[\tag{${\it \%o}_{92}$}\frac{x^5}{5}-\frac{x^4}{4}+\frac{x^3}{3}-\frac{x^2}{2}+x\]

式の整理・簡単化

Maxima では,多項式やいろいろな関数を含んだ式の変形や整理・簡単化が可能です。

展開・因数分解・簡単化

$(x + 1)^2$ を展開します。

In [70]:
(x + 1)^2;
expand(%);
Out[70]:
\[\tag{${\it \%o}_{93}$}\left(x+1\right)^2\]
Out[70]:
\[\tag{${\it \%o}_{94}$}x^2+2\,x+1\]

$x^3 - 1$ を因数分解します。

In [71]:
factor(x^3 - 1);
Out[71]:
\[\tag{${\it \%o}_{95}$}\left(x-1\right)\,\left(x^2+x+1\right)\]

$\displaystyle \frac{3}{x^3 -1}$ を部分分数に展開する例。

In [72]:
partfrac(3/(x^3 - 1), x);
Out[72]:
\[\tag{${\it \%o}_{96}$}\frac{-x-2}{x^2+x+1}+\frac{1}{x-1}\]

ratsimp を使って,上式を「簡単化」します。

In [73]:
ratsimp(%);
Out[73]:
\[\tag{${\it \%o}_{97}$}\frac{3}{x^3-1}\]

三角関数を含む式の簡単化

三角関数や双曲線関数を含む式の簡単化には,trigsimp 関数などを使います。

$\cos^2 x - \sin^2 x$ を「簡単化」します。

In [74]:
trigsimp(cos(x)^2 - sin(x)^2);
Out[74]:
\[\tag{${\it \%o}_{98}$}2\,\cos ^2x-1\]

trigreduce は三角関数同士の積をなるべく減らして「簡単化」します。

In [75]:
trigreduce(cos(x)^2 - sin(x)^2);
Out[75]:
\[\tag{${\it \%o}_{99}$}\frac{\cos \left(2\,x\right)+1}{2}+\frac{\cos \left(2\,x\right)}{2}-\frac{1}{2}\]

通分すればもっと簡単になりますね。このへんが Maxima のつれないところです。

In [76]:
trigsimp(%);
Out[76]:
\[\tag{${\it \%o}_{100}$}\cos \left(2\,x\right)\]

trigexpand 関数は加法定理や倍角の公式を使って「展開」します。

In [77]:
trigexpand(cos(3*x));
trigsimp(%);
Out[77]:
\[\tag{${\it \%o}_{101}$}\cos ^3x-3\,\cos x\,\sin ^2x\]
Out[77]:
\[\tag{${\it \%o}_{102}$}4\,\cos ^3x-3\,\cos x\]

練習

(練習)加法定理・三倍角の公式

  1. 三角関数及び双曲線関数の加法定理を確認しなさい。
  2. $\sin 3x$ を $\sin x$ だけを使って表しなさい。
  3. $\tan 3x$ を $\tan x$ だけを使って表しなさい.
In [78]:
/* 1. のヒント */
/* trigexpand で展開します。*/
trigexpand(sin(x + y))$
trigexpand(cos(x + y))$
trigexpand(tan(x + y))$
trigexpand(sinh(x + y))$
trigexpand(cosh(x + y))$
trigexpand(tanh(x + y))$
In [79]:
/* 2. や 3. のヒント */
trigexpand(sin(3*x));
trigsimp(%);
Out[79]:
\[\tag{${\it \%o}_{109}$}3\,\cos ^2x\,\sin x-\sin ^3x\]
Out[79]:
\[\tag{${\it \%o}_{110}$}\left(4\,\cos ^2x-1\right)\,\sin x\]

$\cos(x)^2$ を $1-\sin(x)^2$ で置き換える例です。

In [80]:
ratsubst(1-sin(x)^2, cos(x)^2, %)$

微分・積分・数値積分

微分・積分を行う関数の表記を以下に示します。

操作 表記例
微分 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)

微分・積分

$x^3 +5x+2$ を $x$ で微分します。

In [81]:
diff(x^3 + 5*x + 2, x);
Out[81]:
\[\tag{${\it \%o}_{112}$}3\,x^2+5\]

$x^3 +5x+2$ を $x$ で 2 階微分します。

In [82]:
diff(x^3 + 5*x + 2, x, 2);
Out[82]:
\[\tag{${\it \%o}_{113}$}6\,x\]

$x^2 \cos x + 2 x \sin x$ を $x$ で積分しま す。

In [83]:
integrate(x^2*cos(x) + 2*x*sin(x), x);
Out[83]:
\[\tag{${\it \%o}_{114}$}2\,\left(\sin x-x\,\cos x\right)+\left(x^2-2\right)\,\sin x+2\,x\,\cos x\]

もう少し簡単になりそうですね。このへんが Maxima のつれないところです。

In [84]:
trigsimp(%);
Out[84]:
\[\tag{${\it \%o}_{115}$}x^2\,\sin x\]

定積分の例:$\displaystyle \int_0^{\pi/2} \sin x \, dx = 1$

以下の例では ' から始めると,その部分の計算を実行せずに形を整えて表示だけしてくれるところを示しています。

In [85]:
'integrate(sin(x), x, 0, %pi/2) = 
 integrate(sin(x), x, 0, %pi/2);
Out[85]:
\[\tag{${\it \%o}_{116}$}\int_{0}^{\frac{\pi}{2}}{\sin x\;dx}=1\]

有名なガウス積分。不定積分はできないのに,無限区間の定積分だけはできる!という例。

In [86]:
'integrate(exp(-x^2), x, -inf, inf) = 
 integrate(exp(-x^2), x, -inf, inf);
Out[86]:
\[\tag{${\it \%o}_{117}$}\int_{-\infty }^{\infty }{e^ {- x^2 }\;dx}=\sqrt{\pi}\]

練習

(練習)関数の傾き

$\displaystyle y = \frac{2x - 1}{x+1}$ が $x \neq -1$ で増加関数であることを示しなさい。

練習

(練習)積分

以下の積分をしなさい。

$$1.\ \int e^x \sin x \ dx \qquad 2. \ \int \frac{x^2 - 2 x - 2}{x^3 -1}\ dx$$

数値積分

解析的に積分できない場合は,数値積分によって近似値を求めることができます。

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

In [87]:
integrate(sin(sin(x)), x, 0, %pi);
Out[87]:
\[\tag{${\it \%o}_{118}$}\int_{0}^{\pi}{\sin \sin x\;dx}\]

解析的に解けないと,そのままで返します。

romberg 関数を使うと,数値的に積分した結果を表示します。

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

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

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

デフォルトの rombergtol の値は...

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

rombergtol を $10^{-12}$ にしてもう一度同じ数値積分を行うと...

In [90]:
rombergtol: 1.e-12$
r2: romberg(sin(sin(x)), x, 0, %pi);
Out[90]:
\[\tag{${\it \%o}_{122}$}1.786487481950052\]
In [91]:
r1 - r2;
Out[91]:
\[\tag{${\it \%o}_{123}$}5.187015705843123 \times 10^{-9}\]

方程式

多項式方程式・連立方程式

Maxima は,式の計算だけでなく,代数方程式を解くことができます.例を以下に示します。

$x^2 +2x−2 = 0$ を $x$ について解き ます.

In [92]:
solve(x^2 + 2*x - 2 = 0, x);
Out[92]:
\[\tag{${\it \%o}_{124}$}\left[ x=-\sqrt{3}-1 , x=\sqrt{3}-1 \right] \]

連立方程式 $x^2 + y^2 = 5, \ x + y = 1$ を $x, y$ について解きます。

In [93]:
solve([x^2+y^2=5, x+y=1], [x, y]);
Out[93]:
\[\tag{${\it \%o}_{125}$}\left[ \left[ x=-1 , y=2 \right] , \left[ x=2 , y=-1 \right] \right] \]

solve の代わりに algsys も使えます。

In [94]:
algsys([x^2+y^2=5, x+y=1], [x, y]);
Out[94]:
\[\tag{${\it \%o}_{126}$}\left[ \left[ x=-1 , y=2 \right] , \left[ x=2 , y=-1 \right] \right] \]

3 次方程式 $x^3 = 8$ を解きます。

In [95]:
solve(x^3 = 8, x);
Out[95]:
\[\tag{${\it \%o}_{127}$}\left[ x=\sqrt{3}\,i-1 , x=-\sqrt{3}\,i-1 , x=2 \right] \]

確認のため,上記の1番目の答えを3乗して簡単化すると...

In [96]:
rhs(%[1])^3, ratsimp;
Out[96]:
\[\tag{${\it \%o}_{128}$}8\]

実数解のみを求める場合は, realroots を使います。

In [97]:
realroots(x^3 = 8);
Out[97]:
\[\tag{${\it \%o}_{129}$}\left[ x=2 \right] \]

練習

(練習)鶴と亀と蟻

鶴と亀と蟻,個体数の合計は 10。足の数は全部で 34 本。蟻は亀より 1 匹少ない。鶴,亀,蟻,それぞれ何羽・何匹?

鶴と亀と蟻の個体数をそれぞれ $x, y, z$ として連立方程式をたてて solve します。蟻の足は何本かわかりますよね?

かつて線形代数の授業でこの問題を出したとき,真っ先に出た質問が「先生,蟻の足は何本ですか?」でした... orz

方程式の数値解

1 変数多項式方程式の数値解を求める関数として, allrootsfind_root があります。

$f(x) := x^3 + 8x - 3$ と定義します。

In [98]:
f(x):= x^3 + 8*x - 3;
Out[98]:
\[\tag{${\it \%o}_{130}$}f\left(x\right):=x^3+8\,x-3\]

allroots を使って,$f(x) = 0$ の全数値解を求めます。

In [99]:
asol: allroots(f(x) = 0);
Out[99]:
\[\tag{${\it \%o}_{131}$}\left[ x=0.3687331875693762 , x=2.846396515370145\,i-0.1843665937846881 , x=-2.846396515370145\,i-0.1843665937846881 \right] \]

realroots を使って $f(x) = 0$ の実 数解を求めます。

In [100]:
rsol: realroots(f(x) = 0);
float(%);
Out[100]:
\[\tag{${\it \%o}_{132}$}\left[ x=\frac{12372633}{33554432} \right] \]
Out[100]:
\[\tag{${\it \%o}_{133}$}\left[ x=0.3687331974506378 \right] \]

allroots で求めた結果と少し違いますね。

精度を指定して再計算してみます。

In [101]:
realroots(f(x) = 0, 1e-16)$
float(%);
Out[101]:
\[\tag{${\it \%o}_{135}$}\left[ x=0.3687331875693762 \right] \]

多項式方程式でない場合には find_root 関数を使って,範囲を指定して解を探します。

関数 $g(x):= (x − 1)\,e^x$ の定義。

$g(x):= (x − 1)\,e^x = 1$ の解は solve などでは求めることができません。

In [102]:
g(x):= (x - 1)*exp(x);
solve(g(x) = 1, x);
Out[102]:
\[\tag{${\it \%o}_{136}$}g\left(x\right):=\left(x-1\right)\,\exp x\]
Out[102]:
\[\tag{${\it \%o}_{137}$}\left[ x=e^ {- x }\,\left(e^{x}+1\right) \right] \]

find_root を使って $0 \leq x \leq 2$ の 範囲で解を探してみます。

出た答えを代入して確かめると...

In [103]:
find_root(g(x) = 1, x, 0, 2);
g(%);
Out[103]:
\[\tag{${\it \%o}_{138}$}1.278464542761074\]
Out[103]:
\[\tag{${\it \%o}_{139}$}1.0\]

練習

(練習)関数の極大値

次の関数が極大値をとる $x$ の値およびその極大値を求めよ。

$$ 1. \ \ f(x) = \frac{x^3}{e^x - 1}, \qquad 2. \ \ g(x) = \frac{x^5}{e^x -1}$$
In [104]:
/* 1. のヒント */
/* f(x) を定義して,それを微分して... */
f(x):= x^3/(exp(x) - 1);
diff(f(x), x);
define(df(x), %);

/* 導関数が 0 になる x を数値的に求める */
find_root(df(x) = 0, x, 1, 5)$

/* 極値は... */
f(%)$
Out[104]:
\[\tag{${\it \%o}_{140}$}f\left(x\right):=\frac{x^3}{\exp x-1}\]
Out[104]:
\[\tag{${\it \%o}_{141}$}\frac{3\,x^2}{e^{x}-1}-\frac{x^3\,e^{x}}{\left(e^{x}-1\right)^2}\]
Out[104]:
\[\tag{${\it \%o}_{142}$}{\it df}\left(x\right):=\frac{3\,x^2}{e^{x}-1}-\frac{x^3\,e^{x}}{\left(e^{x}-1\right)^2}\]

常微分方程式: ode2

Maxima は 1 階および 2 階の常微分方程式も扱うことができます。 まず ode2 関数を使って 解く例を示します。

1 階常微分方程式と初期条件

微分方程式を記述する際は diff の 前の に注意。

ode2 で $\displaystyle \frac{dy}{dt} = a y$ を解きます。

In [105]:
kill(a, t, y)$ /* 念のため,これから使う変数を初期化。 */

ode2('diff(y, t) = a* y, y, t);
Out[105]:
\[\tag{${\it \%o}_{146}$}y={\it \%c}\,e^{a\,t}\]

%c は積分定数です。

1 階常微分方程式の初期条件には ic1 を使い,例えば初期条件を $t = 0$ で $y = 1$ とすると積分定数 %c の値が決まります。

In [106]:
ic1(%, t = 0, y = 1);
Out[106]:
\[\tag{${\it \%o}_{147}$}y=e^{a\,t}\]

2 階常微分方程式と初期条件

2階常微分方程式 $\displaystyle \frac{d^2 x}{dt^2} + K x = 0$ を解く例です。

ここでは単振動の運動方程式を想定しているので $K > 0$です。あとで聞かれないように最初に仮定しておきます。

In [107]:
assume( K > 0);
ode2('diff(x, t, 2) + K*x = 0, x, t);
Out[107]:
\[\tag{${\it \%o}_{148}$}\left[ K>0 \right] \]
Out[107]:
\[\tag{${\it \%o}_{149}$}x={\it \%k}_{1}\,\sin \left(\sqrt{K}\,t\right)+{\it \%k}_{2}\,\cos \left(\sqrt{K}\,t\right)\]

%k1%k2 は積分定数です。

2 階常微分方程式の初期条件は ic2 で設定します。 初期条件を $t = 0$ で $x=x_0, v=v_0$ とすると %k1%k2 が決まります。

In [108]:
ic2(%, t = 0, x = x[0], 'diff(x, t) = v[0]);
Out[108]:
\[\tag{${\it \%o}_{150}$}x=\frac{v_{0}\,\sin \left(\sqrt{K}\,t\right)}{\sqrt{K}}+x_{0}\,\cos \left(\sqrt{K}\,t\right)\]

細かいことですが,x[0] のようにすると Maxima は下付きの添字 $x_0$ として表示します。

常微分方程式: desolve

常微分方程式を解く関数としては desolve もあります。 ode2 と書式が若干違います。

また,初期条件は atvalue 関数で先に設定しておきます。

2 階の常微分方程式 $\displaystyle \frac{d^2 x}{dt^2} + K x = 0$ を解く例を以下に示します。

In [109]:
atvalue(x(t), t = 0, x[0]);
atvalue('diff(x(t), t), t = 0, v[0]);
Out[109]:
\[\tag{${\it \%o}_{151}$}x_{0}\]
Out[109]:
\[\tag{${\it \%o}_{152}$}v_{0}\]
In [110]:
desolve('diff(x(t), t, 2) + K*x(t) = 0, x(t));
Out[110]:
\[\tag{${\it \%o}_{153}$}x\left(t\right)=\frac{v_{0}\,\sin \left(\sqrt{K}\,t\right)}{\sqrt{K}}+x_{0}\,\cos \left(\sqrt{K}\,t\right)\]

ちなみに,$K \rightarrow 0$ の極限で,この解がどうなるかというと...

In [111]:
limit(%, K, 0);
Out[111]:
\[\tag{${\it \%o}_{154}$}x\left(t\right)=v_{0}\,t+x_{0}\]

最初から $K = 0$ とした微分方程式 $\displaystyle \frac{d^2 x}{dt^2} = 0$ を解いても...

In [112]:
desolve('diff(x(t), t, 2) = 0, x(t));
Out[112]:
\[\tag{${\it \%o}_{155}$}x\left(t\right)=v_{0}\,t+x_{0}\]

練習

(練習)重力場中の投げ上げ運動

  1. 高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を求めなさい。運動方程式は, $$ \frac{d^2 y}{dt^2} = -g$$
  2. 速度に比例する空気抵抗がある場合,運動方程式は以下のようになる。 $$ \frac{d^2 y}{dt^2} = -g - \beta \frac{dy}{dt}$$ このとき,高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を 求めなさい。
  3. 前問 2. で求めた解が $\beta \rightarrow 0$ のとき,前問 1. の答えに一致することを示しなさい。

ベクトルと行列

3 次元ベクトルの表記

3 次元ベクトル $\boldsymbol{a} = (a_x,a_y,a_z)$ の表記例。

In [113]:
a: [a[x], a[y], a[z]];
Out[113]:
\[\tag{${\it \%o}_{156}$}\left[ a_{x} , a_{y} , a_{z} \right] \]

$\boldsymbol{a}$ の $x$ 成分(第 1 成分)を取り出すには...

In [114]:
a[1];
Out[114]:
\[\tag{${\it \%o}_{157}$}a_{x}\]

ベクトル $\boldsymbol{b}$ も同様に定義。

In [115]:
b: [b[x], b[y], b[z]];
Out[115]:
\[\tag{${\it \%o}_{158}$}\left[ b_{x} , b_{y} , b_{z} \right] \]

内積

ベクトルの内積 $\boldsymbol{a} \cdot \boldsymbol{b}$ には . (ドット)を使います。

In [116]:
a.b;
Out[116]:
\[\tag{${\it \%o}_{159}$}a_{z}\,b_{z}+a_{y}\,b_{y}+a_{x}\,b_{x}\]

内積を利用して,ベクトルの大きさを返す関数 norm を定義します。

In [117]:
norm(v):= sqrt(v.v);
norm(a);
Out[117]:
\[\tag{${\it \%o}_{160}$}{\it norm}\left(v\right):=\sqrt{v\cdot v}\]
Out[117]:
\[\tag{${\it \%o}_{161}$}\sqrt{a_{z}^2+a_{y}^2+a_{x}^2}\]

外積

ベクトルの外積 $\boldsymbol{a} \times \boldsymbol{b}$ を返す関数 cross を定義します。

In [118]:
cross(u, v):=
           [u[2]*v[3] - u[3]*v[2],
            u[3]*v[1] - u[1]*v[3],
            u[1]*v[2] - u[2]*v[1]]$

成分ごとに見てみます。

In [119]:
cross(a, b)[1];
cross(a, b)[2];
cross(a, b)[3];
Out[119]:
\[\tag{${\it \%o}_{163}$}a_{y}\,b_{z}-b_{y}\,a_{z}\]
Out[119]:
\[\tag{${\it \%o}_{164}$}b_{x}\,a_{z}-a_{x}\,b_{z}\]
Out[119]:
\[\tag{${\it \%o}_{165}$}a_{x}\,b_{y}-b_{x}\,a_{y}\]

練習

(練習)ベクトルの内積,外積

$\displaystyle \boldsymbol{a} = \left(-\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, \sqrt{3}\right), \boldsymbol{b} = \left(-1, 1, -\sqrt{2}\right)$ について,以下の量を計算しなさい。

$$1.\ \ |\boldsymbol{a}|, \quad 2.\ \ |\boldsymbol{b}|, \quad 3. \ \ \boldsymbol{a}\cdot\boldsymbol{b}, \quad 4. \ \ \boldsymbol{a}\times\boldsymbol{b}, \quad 5. \ \ \cos\theta, \quad 6. \ \ \sin\theta$$

ここで $\theta$ は2つのベクトルのなす角。

練習

(練習)ベクトルの 3 重積

  1. スカラー3重積の以下の性質を確かめなさい。 $$\boldsymbol{a}\cdot (\boldsymbol{b}\times\boldsymbol{c}) = \boldsymbol{b}\cdot (\boldsymbol{c}\times\boldsymbol{a}) = \boldsymbol{c}\cdot (\boldsymbol{a}\times\boldsymbol{b}) $$
  2. ベクトル3重積の以下の性質を確かめなさい。 $$\boldsymbol{a}\times (\boldsymbol{b}\times\boldsymbol{c}) = (\boldsymbol{a}\cdot\boldsymbol{c})\,\boldsymbol{b} - (\boldsymbol{a}\cdot\boldsymbol{b})\,\boldsymbol{c} $$
In [120]:
/* 1. のヒント */
kill(c)$

c: [c[x], c[y], c[z]];
a.cross(b, c) - b.cross(c, a)$
expand(%)$
Out[120]:
\[\tag{${\it \%o}_{167}$}\left[ c_{x} , c_{y} , c_{z} \right] \]

行列の簡単な演算

変数を別に再利用する際には,以前の定義を消しておきます。

In [121]:
kill(a, b, c, d)$

2 行 2 列の行列 A の定義。

In [122]:
A: matrix([a, b], [c, d]);
Out[122]:
\[\tag{${\it \%o}_{171}$}\begin{pmatrix}a & b \\ c & d \\ \end{pmatrix}\]

transpose は転置行列。

In [123]:
transpose(A);
Out[123]:
\[\tag{${\it \%o}_{172}$}\begin{pmatrix}a & c \\ b & d \\ \end{pmatrix}\]

determinant は行列式。

In [124]:
determinant(A);
Out[124]:
\[\tag{${\it \%o}_{173}$}a\,d-b\,c\]

逆行列は A^^-1 のように書きます。

In [125]:
A^^-1;
Out[125]:
\[\tag{${\it \%o}_{174}$}\begin{pmatrix}\frac{d}{a\,d-b\,c} & -\frac{b}{a\,d-b\,c} \\ -\frac{c}{a\,d-b\,c} & \frac{a}{a\,d-b\,c} \\ \end{pmatrix}\]

$A^{-1}\, A$ が単位行列になることを確かめます。 行列の積も . (ドット)で 表します。

In [126]:
A^^-1 . A, ratsimp;
Out[126]:
\[\tag{${\it \%o}_{175}$}\begin{pmatrix}1 & 0 \\ 0 & 1 \\ \end{pmatrix}\]

回転行列 O の定義。

In [127]:
O: matrix([cos(x), -sin(x), 0],
          [sin(x),  cos(x), 0],
          [0,      0,       1]);
Out[127]:
\[\tag{${\it \%o}_{176}$}\begin{pmatrix}\cos x & -\sin x & 0 \\ \sin x & \cos x & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}\]

直交行列とは,転置行列が逆行列になっている行列のこと。この回転を表す行列 O が直交行列であることを示してみます。

In [128]:
transpose(O) - O^^-1;
Out[128]:
\[\tag{${\it \%o}_{177}$}\begin{pmatrix}\cos x-\frac{\cos x}{\sin ^2x+\cos ^2x} & \sin x-\frac{\sin x}{\sin ^2x+\cos ^2x} & 0 \\ \frac{\sin x}{\sin ^2x+\cos ^2x}-\sin x & \cos x-\frac{\cos x}{\sin ^2x+\cos ^2x} & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

うーん,この辺が Maxima のつれないところです。

In [129]:
transpose(O) - O^^-1, trigreduce;
Out[129]:
\[\tag{${\it \%o}_{178}$}\begin{pmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

以下のようにしても,直交行列であることが示せます。

In [130]:
O. transpose(O);
Out[130]:
\[\tag{${\it \%o}_{179}$}\begin{pmatrix}\sin ^2x+\cos ^2x & 0 & 0 \\ 0 & \sin ^2x+\cos ^2x & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}\]
In [131]:
trigsimp(%);
Out[131]:
\[\tag{${\it \%o}_{180}$}\begin{pmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}\]

3 次元ベクトル $\boldsymbol{r}$ を,3 行 1 列の行列として定義します。このように成分を縦に並べたものを列ベクトルと言うんでしたね。

In [132]:
r: matrix([x], [y], [z]);
Out[132]:
\[\tag{${\it \%o}_{181}$}\begin{pmatrix}x \\ y \\ z \\ \end{pmatrix}\]

3 次元ベクトルの直交変換の例。

In [133]:
O . r;
Out[133]:
\[\tag{${\it \%o}_{182}$}\begin{pmatrix}x\,\cos x-\sin x\,y \\ \cos x\,y+x\,\sin x \\ z \\ \end{pmatrix}\]

練習

(練習)直交変換されたベクトルの大きさ

$\boldsymbol{r}' \equiv O \,\boldsymbol{r}$ の大きさは $\boldsymbol{r}$ の大きさと同じであることを示しなさい。

In [134]:
/* ヒント */
(O.r) . (O.r)$
trigsimp(%)$

Maxima-Jupyter によるグラフ作成

Maxima-Jupyter を使って,関数のグラフを描くことができます。また,数値データもグラフにすることができます。

In [135]:
/* Jupyter Notebook にインラインでグラフを表示させる場合,*/
/* 最初に以下のようにプロットオプションをセット。 */
set_plot_option([svg_file, "maxplot.svg"])$

2 次元グラフ plot2d の例

2 次元グラフ $y = f(x)$ を描く一般的な表記例は以下の通りです。

plot2d(f(x), [x, a, b])

例として,$y = \sin x$ のグラフを $−2\pi \leq x \leq 2\pi$ の範囲で描きます。基本的な定数の一つである円周率 $\pi$ は Maxima では %pi と書くのでしたね。

In [136]:
plot2d(sin(x), [x, -2*%pi, 2*%pi])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 sin(x) x gnuplot_plot_1 -1 -0.5 0 0.5 1 -6 -4 -2 0 2 4 6

3 次元グラフ plot3d の例

3 次元グラフ $z = g(x, y)$ を描く一般的な表記例は以下の通りです。

plot3d(g(x,y), [x, a, b], [y, c, d])

例として,$z = x \sin y$ のグラフを $-3 \leq x \leq 3, -4 \leq y \leq 4$ の範囲で描きます。

In [137]:
plot3d(x*sin(y), [x, -3, 3], [y, -4, 4])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -3 -2 -1 0 1 2 3 -4 -3 -2 -1 0 1 2 3 4 -3 -2 -1 0 1 2 3 gnuplot_plot_1 x*sin(y) x y z

数値データ・リストのグラフ作成例

Maxima は,リスト形式のリスト形式の数値データもグラフ化することができます。

以下の例では,最初(左側)の数値を $x$ 座標の値,次(右側)の数値を $y$ 座標の値として 6 個の点 からなるリスト data の各点をつないだ折れ線グラフを描きます。

In [138]:
data: makelist([i, i^2], i, 0, 5);
plot2d([discrete, data])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x gnuplot_plot_1 0 5 10 15 20 25 0 1 2 3 4 5
Out[138]:
\[\tag{${\it \%o}_{188}$}\left[ \left[ 0 , 0 \right] , \left[ 1 , 1 \right] , \left[ 2 , 4 \right] , \left[ 3 , 9 \right] , \left[ 4 , 16 \right] , \left[ 5 , 25 \right] \right] \]

オプションを設定してプロットする例。

ここでは,[style, points]を設定して線でつながずに点で描きます。プロットオプションの style には,points の他にも,lineslinespointsdots が設定できます。

In [139]:
plot2d([discrete, data], [style, points])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x gnuplot_plot_1 0 5 10 15 20 25 0 1 2 3 4 5

でかい点ですね。点のサイズを変更して描きます。

In [140]:
plot2d([discrete, data], [style, [points, 1.5]])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x gnuplot_plot_1 0 5 10 15 20 25 0 1 2 3 4 5

複数のグラフを重ねて表示

Maxima で複数のグラフを重ねて表示する例を示します。

$y = x^2 - 1$ と $y = 4 x - 5$ の 2 つのグラフを重ねて描きます。

$x$ の範囲は $-5 \leq x \leq 5$ です。 以下の例でわかるように,重ねて描きたい関数を [x^2 - 1, 4*x - 5] のように [] で囲んだリスト形式にして plot2d コマンドを使います。

In [141]:
plot2d([x^2 - 1, 4*x - 5], [x, -5, 5])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 x x^2-1 x^2-1 4*x-5 4*x-5 -25 -20 -15 -10 -5 0 5 10 15 20 25 -4 -2 0 2 4

いくつかオプションを設定して描く例です。

以下の例では,縦軸 ($y$ 軸) の表示範囲を制限し たために,

plot2d: some values were clipped. というメッセージが出ています。

In [142]:
plot2d([x^2 - 1, 4*x - 5], [x, -5, 5],
       [y, -5, 10],                         /* 縦軸の表示範囲の設定 */
       [color, red, black],                 /* 線の色の設定 */ 
       [style, [lines, 4], lines],          /* 片方の線の太さを4に */
       [gnuplot_preamble, "set key bottom"] /* 凡例の位置 */
       )$
plot2d: some values were clipped.
plot2d: some values were clipped.
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 x x^2-1 x^2-1 4*x-5 4*x-5 -4 -2 0 2 4 6 8 10 -4 -2 0 2 4

上のグラフをみると,2 つのグラフは 1 点で接しているように見えます。

実際に $y = 4x - 5$ が接線であることを示してみましょう。

$x^2 - 1 = 4 x - 5$ の解が1つであることを示せば十分ですね。

In [143]:
sol: solve(x^2 - 1 = 4*x - 5, x);
Out[143]:
\[\tag{${\it \%o}_{194}$}\left[ x=2 \right] \]

$y = f(x)$ の $x = x_1$ における接線の方程式は, $$ y - f(x_1) = f'(x_1) \,(x - x_1)$$ ですから...

In [144]:
f(x):= x^2 - 1;
define(df(x), diff(f(x), x));

x1: rhs(sol[1]);
y = df(x1) * (x - x1) + f(x1), ratsimp;
Out[144]:
\[\tag{${\it \%o}_{195}$}f\left(x\right):=x^2-1\]
Out[144]:
\[\tag{${\it \%o}_{196}$}{\it df}\left(x\right):=2\,x\]
Out[144]:
\[\tag{${\it \%o}_{197}$}2\]
Out[144]:
\[\tag{${\it \%o}_{198}$}y=4\,x-5\]

できたら,$y = x2 - 1$ の $x = 1$ 等における接線のグラフも描いてみましょう。

数値データファイルを読み込む

Maxima では,あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。

以下のような内容のファイル test.dat がカレントディレクトリにあるとします。

0   0
1   1
2   4
3   9
4   16
5   25

データファイルを読み込む Maxima の関数は,read_nested_list を使います。

In [145]:
dat: read_nested_list("test.dat");
Out[145]:
\[\tag{${\it \%o}_{199}$}\left[ \left[ 0 , 0 \right] , \left[ 1 , 1 \right] , \left[ 2 , 4 \right] , \left[ 3 , 9 \right] , \left[ 4 , 16 \right] , \left[ 5 , 25 \right] \right] \]
In [146]:
plot2d([discrete, dat],
       [x, 0, 6], [y, 0, 28],                   /* x, y の表示範囲を設定 */ 
       [style, [points, 1]], [color, red],      /* 大きさ,色 */
       [gnuplot_preamble, "set key top left"],  /* 凡例の位置 */ 
       [legend, "データ"]                       /* 凡例表示の設定 */
       )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x データ データ 0 5 10 15 20 25 0 1 2 3 4 5 6

数値データと理論曲線を重ねて表示

前節の数値データファイル test.dat と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。

In [147]:
plot2d([[discrete, dat], x^2], [x, 0, 5], [y, 0, 28],
       [style, [points, 1.5], lines], [color, red, blue],
       [gnuplot_preamble, "set key top left"], 
       [legend, "データ", "理論曲線 y=x^2"]
       )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 x データ データ 理論曲線 y=x^2 理論曲線 y=x^2 0 5 10 15 20 25 0 1 2 3 4 5

媒介変数表示の 2 次元グラフ

半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です.この円を描くには,$y = \pm \sqrt{1 - x^2}$ として $y = f(x)$ の形にするよりも,以下のような媒介変数表示にしたほうが簡単でしょう。

$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$

媒介変数表示の 2 次元グラフを Maxima で描くには以下のようにします。

In [148]:
plot2d([parametric, cos(t), sin(t), [t, 0, 2*%pi]])$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 sin(t) cos(t) gnuplot_plot_1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1

グラフの縦横の 比を1:1にして円らしく見えるように,[same_xy] を設定してみます。

[yx_ratio, 1] でも同等です。

In [149]:
plot2d([parametric, cos(t), sin(t), [t, 0, 2*%pi]],
        [same_xy],
        [xlabel, "x"], [ylabel, "y"], [legend, "x^2+y^2=1"],
        [x,-1.1, 1.1], [y, -1.1, 1.1]
      )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x x^2+y^2=1 x^2+y^2=1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1

練習

(練習)楕円のグラフ(楕円中心を原点にして)

同様にして,楕円のグラフを描くこともできます。 原点を中心とし,長半径 $a$,単半径 $b$ の楕円の式は

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

媒介変数表示では,$x=a \cos t, \ y=b \sin t, \ (0\leq t\leq 2\pi)$と書けます。 $a$,$b$ に適当な値を入れて楕円のグラフを描きなさい。

In [150]:
/* ヒント */
x(a, t):= a*cos(t)$
y(b, t):= b*sin(t)$

plot2d([parametric, x(1,t), y(1,t), [t, 0, 2*%pi]],
        [same_xy],
        [xlabel, "x"], [ylabel, "y"], [legend, "a=1, b=1"],
        [x,-2, 2], [y, -2, 2]
      )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x a=1, b=1 a=1, b=1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

海王星と冥王星の軌道

焦点を原点とした楕円のグラフ

太陽からの万有引力を受けて運動する惑星は,太陽(二体問題では質量中心)を焦点の一つとした楕円運動を描いて運動します。 焦点の一つを原点とし,長半径 $a$,離心率 $e$ の楕円の方程式は,極座標 $r$,$\varphi$ を使って以下のように表すことができます。

$$ r= \frac{a (1−e^2)}{ 1 + e\cos\varphi}$$

さて,かつては第9惑星,現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $𝑎_P=39.767\mbox{AU}$,離心率 $𝑒_P=0.254$ を使って冥王星の軌道を描きます。

まず,極座標表示の楕円の式を関数として定義します。

In [151]:
r(a, e, t):= a*(1-e^2)/(1+e*cos(t))$
x(a, e, t):= r(a, e, t)*cos(t)$
y(a, e, t):= r(a, e, t)*sin(t)$
In [152]:
aP: 39.767$
eP: 0.254$
plot2d([parametric, x(aP, eP, t), y(aP, eP, t), [t, 0, 2*%pi]],
        [same_xy],
        [xlabel, "x"], [ylabel, "y"], [legend, "冥王星"],
        [x,-50, 50], [y, -50, 50]
      )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x 冥王星 冥王星 -40 -20 0 20 40 -40 -20 0 20 40

では,この冥王星の軌道と,その内側をまわる海王星の軌道を重ねて描いてみましょう。

海王星の軌道長半径は $a_N=30.1104\mbox{AU}$,離心率は $e_N=0.0090$ と小さいので簡単のために $e_N=0$ として扱います。実際の2つの天体の軌道は同じ平面上にありませんが,太陽からの距離のみをみるために,ここでは同一平面上に描きます。

In [153]:
aN: 30.1104$
eN: 0$
plot2d([[parametric, x(aP, eP, t), y(aP, eP, t), [t, 0, 2*%pi]],
        [parametric, x(aN, eN, t), y(aN, eN, t), [t, 0, 2*%pi]]
       ],
        [same_xy],
        [xlabel, "x"], [ylabel, "y"], [legend, "冥王星", "海王星"],
        [x,-50, 50], [y, -50, 50]
      )$
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 y x 冥王星 冥王星