大学ではMathimaticaなどで方程式を解くことができた。
Mathematicaのフリー版、とも言える"maxima"というソフトを知ったので、いろいろ試した。
maximaのプロンプトから操作する方法はネットで調べることができたが、シェルスクリプトから触る方法はヒットしなかった。
いわゆる"数値計算"で系をシミュレートしたい時には、面倒である。
(maximaでfor文が使えるらしいことを見つけた。が、覚えるのが・・・）
なのでここでは、シェルスクリプトを使って、maximaに方程式を解かせる方法について、忘備録を残す。
maxima自体の使いかたは、他の方のページを参考にする。
Maximaで方程式・連立方程式を解く、数列を求める
前提
maximaがインストールしてあること。
Macでは、
brew install maxima
で入った。
シェルから操作する限りは不要だが、maximaのGUIである、wxmaximaも入れると良いと楽しいと思う。
maxima ファイルの作り方
$ a x + b = 0 $
を$x$ について解くのであれば、maximaの入力ファイルはsolve()を使って、
solve(a*x+b=0,x);
とする。
例えば、$ 2 x - 3 y = 0 $ であれば、
a:2;
b:-3;
solve(a*x+b=0,x);
とすれば良い。
maxima ファイルの読み込み方
maxima -b file.max
maxima -q -b file.max
maxima -q -batch=file.max
maxima --very-quiet -batch=file.max
-b or --batch=でファイル名の指定
-qや--very-quietとすると、いらない出力が減っていく
シェルスクリプトからmaximaを触る
例1：一次方程式の解
set a = $argv[1]
set b = $argv[2]
set filename=$$.max
cat << EOF > $filename
a:$a;
b:$b;
solve(a*x+b=0,x);
EOF
maxima --very-quiet --batch=$filename
rm $filename
ヒアドキュメントcat << EOF > $filename〜EOFを使って、maxima入力ファイルを作っている。
こうすることで、シェルスクリプトから変数を渡すことができる。
最終行では、ヒアドキュメントで作ったmaxima入力ファイルを削除している。
(シェルを動かした後に、何も残さない様に・・・）
上記のスクリプトの使い方は、
例えば、$ 2 x - 3 y = 0 $ であれば、
csh solve.csh 2 -3
とすれば良い。
解は
solve(a*x+b = 0,x)
3
[x = -]
2
と出力してくれる。
もしも、浮動小数点表示したいなら、float()関数にsolve()式を入れれば良い。
すなわち、
float(solve(a*x+b=0,x));
とする。
例2：円と直線の交点
円周 $x^2+y^2=c^2$ と直線 $y=ax+b$ の交点を求めるスクリプト
set a = $argv[1]
set b = $argv[2]
set c = $argv[3]
set filename=$$.max
cat << EOF > $filename
a:$a; b:$b; c:$c;
solve( [x^2+y^2=c^2, y=a*x+b], [x,y] );
EOF
maxima --very-quiet --batch=$filename
例えば、直線$y=2x+1$ と $x^2+y^2=1$の交点は、
csh solve2.csh 2 1 1
とすれば、
solve([x^2+y^2 = c^2,y = a*x+b],[x,y])
4 3
[[x = - -, y = - -], [x = 0, y = 1]]
5 5
と、教えてくれる。
# 高校生なら、これで課題も楽チンですね。
ところで、"解なし”になるはずの直線$y=2x+5$ と $x^2+y^2=1$の交点は、
csh solve2.csh 2 5 1
とすれば、
solve([x^2+y^2 = c^2,y = a*x+b],[x,y])
2 sqrt(5) %i + 10 4 sqrt(5) %i - 5
[[x = - -----------------, y = - ----------------],
5 5
2 sqrt(5) %i - 10 4 sqrt(5) %i + 5
[x = -----------------, y = ----------------]]
5 5
と、出力される。
%iは虚数$i$を表していて、虚数解を教えてくれる。
Tips
エスケープシーケンス
load("mnewton")\$
load("mnewton")$など、$を使うときは、エスケープシーケンス\を使うこと。
maximaの出力の整形
(%o12) [[A = 1.17e-8, B = 3.98e-9]]
これは、二次元配列(maximaではリストというらしい）に格納されている。
いったん、
ans:%[12]
など、リストansに格納してから取り出すと吉。
ans[1][1]には "A = 1.17e-8"、
ans[1][2]には、"B = 3.98e-9"が入っている。
rhs()関数は、右辺のみを取り出す関数なので、rhs(ans[1][1])とすれば、1.17e-8だけが取り出せる。
つまり、結果だけを表示したければ、
print(rhs(ans[1][1]),rhs(ans[1][2]) ;
とすれば良い。