第2版 update: 2021年10月 葛西 真寿 弘前大学大学院理工学研究科
Maxima-Jupyter を使った Maxima プログラミング入門。 はじめての Python プログラミングの内容を Maxima で書いてみたものです。
Maxima はコンピュータ代数システム(いわゆる数式処理システム)ですが,構造化プログラミングの機能もあり,条件分岐や繰り返し処理などを利用することで,プログラミング言語としても利用できます。
エンゲル係数とは,家計の消費支出に占める飲食費(食料費,食費)の割合(パーセント単位)です。
消費支出 110,000円,飲食費 37,000円の人のエンゲル係数を求める例。Maxima では以下のような電卓的な使い方ですぐ答えが出ます。
本稿のように Jupyter Notebook 環境では,以下の1行を書いたあとに,Shift キーを押しながら Enter キー(または Return キー)を押して実行します。
Maxima の文末は,;
か $
で終わります。$
で終わる場合は結果が表示されません。
37000/110000 * 100;
Maxima では割り算の結果は厳密な有理数として出力されます。(近似的な)実数(小数点数)で表示させたい場合は,とりあえずは数値のどれか1つにでも小数点をつけます。(あとで float()
関数で実数に変換することも説明します。)
37000/110000 * 100.0;
ここからは,単に電卓的に使うだけではなく,応用・発展も見据えて以下のようなプログラム文を作成して実行することにします。
以下の例では,
/*
と */
で囲まれる部分はコメント。コメントはそれを読む人間のための注釈であり,プログラムの実行時には無視されます。spending
に消費支出である 110000
(円)という数値を代入し,food
に飲食費である 37000
(円)という数値を代入し,eng
に代入し,表示させています。Maxima では,変数への代入は :
です。
/* エンゲル係数を求める
消費支出 spending, 飲食費 food, エンゲル係数 eng */
spending: 110000$
food: 37000$
eng: food / spending * 100.0;
意に反して,プログラムがうまく動作しない場合があります。たとえば,上記の 3行目を以下のように書いた場合。
掛け算を示す演算子は *
,割り算を示す演算子は /
であるところを誤って,x
や \
と書くと,これらの記号は Maxima では演算子として定義されていないのでエラーとなります。このようなエラーは文法エラー(incorrect syntax) と呼ばれ,Jupyter Notebook では,以下のように誤りやその箇所を指摘してくれます。
spending: 110000$
food: 37000$
eng: food x 100 \ spending;
別の種類の誤りもあります。たとえば,パーセントにするためにかける「100.0」を誤って「1000.0」と入力した場合...
spending: 110000$
food: 37000$
eng: food / spending * 1000.0;
上記の例では incorrect syntax
(文法エラー)の警告が出ませんが,答えは明らかに変です。(エンゲル係数の最大値は100だから。)
このような場合はコンピュータにとって誤りではなく,指示通りに計算し,その結果を表示します。 このようなプログラム作成者の不注意や勘違いによって生じるエラーには注意が必要です。
プログラムが正しく動作しない場合は,誤りの箇所を探し,修正して再度動作を確認します。正しく動作するようになるまで,プログラムを修正し,実行(Shift + Enter)します。この作業をデバグといいます。
上記の例の 1行目では,spending という名前の変数に 110000 という 数値を代入しています。変数名はアルファベットで始まる英数字列にします。大文字と小文字は区別されます。
データ型には数値(整数・有理数・浮動小数点数)と文字列などがあります。
変数 moji に文字列を代入する場合は,以下のように " で囲みます。
moji: "エンゲル係数"$
print(moji)$
2 + 50 - 5 * 6 /2;
2 + (50 - 5 * 6)/2;
文字列の連結演算は concat()
を使います。以下の例を参照。
mo : "エンゲル"$
ji : "係数"$
moji : concat(mo, ji);
文字列と数値を直接連結するときも,concat()
を使います。
spending: 110000$
food: 37000$
eng: food / spending * 100.0$
output : concat("エンゲル係数は ", eng, " です。")$
print(output)$
文字列と連結しなくても,以下のようにも書くことができます。
print()
文は,カンマ (,
) で区切ることによって複数の項目(変数)を表示することができます。
print("エンゲル係数は ", eng, " です。")$
Maxima の割り算は以下のように厳密な値を返します。
17 / 3;
浮動小数点数表示する場合は,すでに述べたように数値に小数点をつけて整数でないことをアピールするか,浮動小数点数に変換する関数 float()
を使います。
17.0 / 3;
float(17 / 3);
べき乗は **
(または ^
)です。以下の例では $2^4 = 2 \times 2 \times 2 \times 2$ と $2^{1/2} = \sqrt{2}$ の値を表示します。
2**4;
2**(1 / 2);
float(2**(1 / 2));
1 天文単位(1 au)の距離を光速 c で進むのに要する時間はいくらか。ただし,
149597870700
m299792458
m/s である。参考:
上記の例では,プログラムの中で消費支出や飲食費を決めていましたが,今度はプログラムを実行する人が入力できるようにしてみます。
read()
関数は入力された値を返します。
注意: read()
に対して値をキーボードから入力する場合も,最後は ;
をつけてから Enter キーを押します。
spending : read("消費支出? ")$
food : read("飲食費? ")$
eng : float(food / spending * 100)$ /* float() で変換 */
print("エンゲル係数は ", eng, " です。")$
(これまでの例では1回しか計算していませんが)よく使う計算を「関数」として定義する例です。
以下の例では,1行目でエンゲル係数を計算する関数 calc()
を定義し,5行目でその関数を呼び出しています。
read()
関数に足してキーボードから数値を入力したら最後に ;
をつけるのを忘れないで!
calc(spending, food):= float(food / spending * 100)$
spending : read("消費支出? ")$
food : read("飲食費? ")$
eng : calc(spending, food)$
print("エンゲル係数は ", eng, " です。")$
いったん定義された関数は,以下のように何回でも利用できます。
calc(120000, 53000);
消費支出を入力すると,エンゲル係数が 20 ~ 40 となる飲食費を表示するプログラムをつくりなさい。
Maxima には,$\sin x, \cos x, \tan x$などの数学関数が組み込まれています。また,数学定数としては円周率 %pi
などが使えます。
pi: %pi$
print(pi, sin(pi/2), cos(pi/3), tan(pi/4))$
以下の例では,x
以下の最大の整数を与える関数 floor()
を使って,小数点以下を切り捨てたエンゲル係数の値を表示させています。
また,エンゲル係数を計算する部分はすでに関数 calc()
として定義していますが,プログラム全体を関数 engel()
として定義してしまいます。
以下の例のように,block()
は複数の関数や式を ,
(カンマ)で区切って記述します。
engel():= block(
spending : read("消費支出? "),
food : read("飲食費? "),
eng : calc(spending, food),
print("エンゲル係数は ", floor(eng),
" です。(小数点以下切り捨て)")
)$
engel()$
Maxima で使える数学関数や数学定数については,以下のマニュアルのページを参照してください。
小数点以下を切り上げたエンゲル係数の値を表示させるように上記のプログラムを変更しなさい。
ヒント: Maxima 5.42.2 Manual: 10. Mathematical Functions の ceiling()
を参照。
Maxima の組み込み関数に round()
があります。数値を丸める関数ですが,いわゆる四捨五入とはちょっと違うようです。
以下の例からわかるように,小数点以下が 0.5
未満の場合と 0.5
を超える場合は四捨五入と同様の丸め方ですが,ちょうど 0.5
のときは,偶数に丸められていることがわかります。
これを偶数への丸めといい,端数が 0.5
より小さいなら切り捨て,端数が0.5
より大きいなら切り上げ,端数がちょうど 0.5
なら切り捨てと切り上げのうち結果が偶数となる方へ丸めます。
print(round(12.49), " ", round(12.50), " ", round(12.51))$
print(round(13.49), " ", round(13.50), " ", round(13.51))$
engel():= block(
spending : read("消費支出? "),
food : read("飲食費? "),
eng : calc(spending, food),
if eng >= 40 then
print("エンゲル係数は ", round(eng),
" です。高い値です。")
else
print("エンゲル係数は ", round(eng), " です。")
)$
engel()$
engel()$
if
文の書式は以下の通りです。
if 条件式1 then
条件式1が満たされる場合に実行する文
elseif 条件式2 then
条件式2が満たされる場合に実行する文
else
上記以外の場合に実行する文 $
上の例では,eng >= 40
つまり eng
が 40
以上のとき,という関係演算子を使いました。他の例は以下のとおりです。
a < b
a
が b
より小さい
a <= b
a
が b
以下
a > b
a
が b
より大きい
a >= b
a
が b
以上
is(equal(a, b))
は a
と b
が等しいとき true, そうでないとき false を返します。equal(a, b)
の否定は notequal(a, b)
です。
以下の例では,1行目で代入された消費支出に対してエンゲル係数を計算します。
1回だけ計算して表示するのではなく,2行目の飲食費の初期値(ここでは 40000
円)から,5行目の条件式(ここでは 80000
円以下)が成り立つ場合に,8行目にあるように 5000
円ずつ増やしながらエンゲル係数を表示します。
spending : 110000 $
food : 30000 $
print("消費支出", spending, "円に対するエンゲル係数の値")$
while food <= 50000 do(
eng : calc(spending, food),
print(" 飲食費", food, "円の場合: ", eng),
food : food + 5000
)$
while
文の書式は以下の通りです。
while 条件式 do(
条件式が成り立つ場合に実行する式1,
条件式が成り立つ場合に実行する式2, ... ,
条件式が成り立つ場合に実行する式
)
同様の繰り返し処理を for
文を使って行うこともできます。
spending : 110000 $
food : 30000 $
print("消費支出", spending, "円に対するエンゲル係数の値")$
for i: 30000 thru 50000 step 5000
do(
eng : calc(spending, food),
print(" 飲食費", food, "円の場合: ", eng),
food : food + 5000
)$
for
文の書式は以下の通りです。
for ivar: istart thru iend step istep
do(
ivar が istart から iend の間,実行される文
)$
ivar
は整数変数で,その初期値は istart
。文が実行されるごとに,ivar
の値は istep
だけ加算され,iend
の値になるとループが終了します。
step
の部分を省略すると,1
ずつ加算されます。
以下の値を求めるプログラムを作成せよ。
$$\sum_{n = 1}^{10} n^2$$a : [5, 4, 3, 2, 1]$
print("全成分を一挙に表示")$
print(a)$
少し表示を工夫してみます。sprint()
関数は(改行せずに)横に続けて出力します。
print("全成分を横に表示")$
for i thru length(a) do sprint(a[i])$
print("全成分を縦に表示")$
for i thru length(a) do print(a[i])$
5個の成分を全てゼロにしてリストを初期化する例です。
/* リスト b の5つの要素を全て 0 に*/
b : makelist(0, 5)$
print(b)$
/* 3番目の要素の値を設定 */
b[3] : 3$
print(b)$
以下のようなテキストファイル mydata.txt
を用意します。
$ cat mydata.txt
1
2
3
4
5
6
ファイル mydata.txt
の数値をリスト a
に読み込みます。
a : read_list("mydata.txt");
/* 要素の総和,平均,最大値,最小値を出力します。*/
print("要素の総和は", lsum(i, i, a))$
print("要素の平均は", float(lsum(i, i, a)/length(a)))$
print("最大値は ", lmax(a))$
print("最小値は ", lmin(a))$
次は,以下のような「6行2列」のデータが入ったテキストファイル mydata2.txt
を読み込む例です。
$ cat mydata2.txt
1 1
2 4
3 9
4 16
5 25
6 36
b : read_nested_list("mydata2.txt")$
print(b)$
/* 入れ子になっているリストの成分を取り出す方法 */
b[2];
b[2][2];
リストではなく,行列として読み込む場合は以下のようにします。
M : read_matrix("mydata2.txt");
/* 行列として読み込んだ場合の成分の取り出し方 */
M[2, 2];
ファイルへの書き込み例です。この例では,$𝑥$ と $\sin 𝑥$ の値を $𝑥=0.1$ から $𝑥=1.0$ まで書き込んでいます。
x : makelist([0.1 * i, sin(0.1 * i)], i, 1, 10);
write_data(x, "myoutput.txt")$
$ cat myoutput.txt
0.1 0.09983341664682815
0.2 0.1986693307950612
0.3 0.2955202066613396
0.4 0.3894183423086505
0.5 0.479425538604203
0.6000000000000001 0.5646424733950355
0.7000000000000001 0.6442176872376911
0.8 0.7173560908995228
0.9 0.7833269096274834
1.0 0.8414709848078965
行列からリストへの変換は args()
関数で行います。
M;
matrixp(M);
m : args(M);
listp(m);
/* リストと行列で成分の参照法が異なる */
M[6,2];
m[6][2];
また,リストから行列への変換は以下のようにします。
b;
listp(b);
M1: apply('matrix, b);
matrixp(M1);
/* リストと行列で成分の参照法が異なる */
b[6][2];
M1[6,2];
実数を表示する際,全体の文字数や小数点以下の桁数を指定して表示させる例です。
例えば "~8,5f"
は全体で 8 桁の表示範囲をとり,小数点以下は 5 桁表示します。"~%"
は改行を表します。
sqrtwo : float(sqrt(2))$
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~f~%", sqrtwo)$
printf(true, "~8,5f~%", sqrtwo)$
printf(true, "~18,15f~%", sqrtwo)$
printf(true, "2の平方根は ~18,15f です~%", sqrtwo)$
全体の桁数のみを指定することもできます。例えば,"~7f"
のように全体の桁数 7 桁のみを指定します。また,小数点以下の桁数のみを"~,2f"
のように指定することもできます。
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~7f~%", sqrtwo)$
printf(true, "~,2f~%", sqrtwo)$
以下の出力を,エンゲル係数の値を(四捨五入して)小数点以下1桁まで表示するように変更しなさい。
spending : 110000$
food : 37000$
eng : float(food * 100 / spending)$
print("エンゲル係数は ", eng, " です。")$
非常に大きい数や小さい数を表示するときには, $1.9891 \times 10^{30}$ や $6.67430 \times 10^{-11}$ などのような指数表示をします。
例えば "~9,2e"
では,全体で 9 桁の表示分を確保し,小数点以下は 2 桁表示します。
c : 299792458$
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~9,2e~%", c)$
printf(true, "~15,8e~%", c)$
printf(true, "光速は ~15,8e m/s です。~%", c)$
全体の桁数のみを指定したり,小数点以下の桁数のみを指定することもできます。
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~,2e~%", c)$
printf(true, "~15e~%", c)$
文字列の書式指定は "~a"
です。また,以下では substring()
で文字列の一部を取り出す処理を行なう例を示します。
substring(str, i, j)
は文字列 str
の i
番目から始まり,j
番目以降を切り取った(つまり,i
番目から j-1
番目までの)文字列を与えます。
ten : "1234567890"$
printf(true, "~a~a 桁数~%", ten, ten)$
printf(true, "~a~a~%", substring(ten, 1, 6),
substring(ten, 1, 6))$
printf(true, "~a~a~%", substring(ten, 1, 7),
substring(ten, 1, 5))$
整数の書式指定は "~d"
です。"~5d"
のように桁数を指定することもできます。
v : [1, 23, 456, 7890]$
for i thru length(v) do print(v[i])$
/* 書式指定で表示桁数を5桁とり,右寄せで整数を表示 */
for i thru length(v) do printf(true, "~5d~%", v[i])$
m: [[1, 0, 0, 0],
[0, 23, 0, 0],
[0, 0, 456, 0],
[0, 0, 0, 7890]];
/* sprint() は改行せずに出力。newline() は改行のみ。*/
for i thru length(m)
do(for j thru length(m[i]) do sprint(m[i][j]),
newline()
)$
/* 書式指定で表示桁数を5桁とり,右寄せで各成分を表示 */
for i thru length(m)
do(for j thru length(m[i]) do printf(true, "~5d" ,m[i][j]),
newline()
)$
この Notebook では,プログラミング言語として共通の部分に焦点を絞った Maxima の使い方について解説しました。
微分・積分や方程式の解法,さらにはグラフ作成といった,コンピュータ代数システム Maxima 固有の機能の活用例については,以下のページを参考にしてください。
また,Maxima のマニュアルは