シェルスクリプトからmaximaで方程式を解かせる

大学ではMathimaticaなどで方程式を解くことができた。
Mathematicaのフリー版、とも言える"maxima"というソフトを知ったので、いろいろ試した。

maximaのプロンプトから操作する方法はネットで調べることができたが、シェルスクリプトから触る方法はヒットしなかった。
いわゆる"数値計算"で系をシミュレートしたい時には、面倒である。
(maximaでfor文が使えるらしいことを見つけた。が、覚えるのが・・・）


なのでここでは、シェルスクリプトを使って、maximaに方程式を解かせる方法について、忘備録を残す。


maxima自体の使いかたは、他の方のページを参考にする。
[Maximaで方程式・連立方程式を解く、数列を求める](https://pianofisica.hatenablog.com/entry/2018/12/04/044834)

# 前提

maximaがインストールしてあること。
Macでは、

```shell
brew install maxima
```

で入った。
シェルから操作する限りは不要だが、maximaのGUIである、`wxmaxima`も入れると良いと楽しいと思う。

# maxima ファイルの作り方

$ a x + b = 0 $

を$x$ について解くのであれば、maximaの入力ファイルは`solve()`を使って、

```file.max
solve(a*x+b=0,x);
```

とする。

例えば、$ 2 x - 3 y = 0 $　であれば、

```file.max
a:2;
b:-3;
solve(a*x+b=0,x);
```

とすれば良い。


# maxima ファイルの読み込み方

```shell
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：一次方程式の解

```shell:solve.csh
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$ の交点を求めるスクリプト

```shell:solve2.csh

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]) ;
```
とすれば良い。




# 参考


[Maximaで方程式・連立方程式を解く、数列を求める](https://pianofisica.hatenablog.com/entry/2018/12/04/044834)