数式処理ソフトのすすめ

挨拶
お前は誰だ？
真面目に
数式処理ソフトとは？
Maximaを導入する
課題は楽に解ける！！！
LaTeX出力について
終わりに

挨拶

こんにちは。bokuroroです。この記事はUEC Advent Caleder2021の20日目の記事です。

この間初めてはてなブログのアカウントを作ったのですが、文字を書くって大変ですね(クソ雑魚)

あまりブログ文化に親しんで来なかったので練習だと思ってやってみます。いまだに文字を書くのって苦手なんですよね...(遠い目)

ところで昨日はsuzukeくんの「Wiki、始めました」という記事でした。個人用で使えるwiki作って情報を整理したいなと思っていたのですが、全然知識がないのと時間がないのとで、なんやかんややってないですね。というかVLLとっくの昔にwikiを導入してバリバリ情報共有していると思ったので、最近wikiが出来たのが少し意外でしたね。引き継ぎ資料とかwiki上に全部まとめたらますますの発展が期待できそうです。

suzuke.dev

お前は誰だ？

僕が誰なのかはTwitterをみればわかります？

twitter.com

知らない人もいるかもなので自己紹介しておくと、III類の電子工学プログラムに所属している3年生です。MMAに所属しています。Advent Calender、今日からMMAメンツが連続して並んでいて面白いですね。Twitterを見ればわかると思いますが、にわかオタクでとあるとジョジョが好きです。推しはとある科学の超電磁砲の黒子が好きですね。黒子が一番可愛いんですよ。

最近研究室配属決まりましたわーい、叫んでもいいですか？　グギュグバァッ！！！(ディアルガの鳴き声)

youtu.be

(ちなみにダイパリメイクは買う予定ない)

真面目に

この記事は数式処理ソフトを導入して楽しい課題消化ライフを送ろう！！というものです。皆さんは、数式をまだ手計算でやっているのですか？ 手計算でやると、線形代数の固有値、固有ベクトルを求めたりとか、微分方程式を解くのは死ぬほど辛いですよね... 計算量も多いし、その分計算ミスも多くなってきますし... そこで数式処理ソフトを導入して少しでも課題を楽にしようという趣旨の記事です。

数式処理ソフトとは？

C言語やpythonといったプログラミング言語は数値や四則演算を扱うのは得意ですが、xやyといった式を処理したり、積分したり、方程式を扱うのは難しいです。そこで、それを自動で行なってくれるのが数式処理ソフトです。数式処理ソフトの例としては、mathematicaやMaple,Maximaといったものが挙げられます。大学で無料で使えるのはMapleですね。

www.cc.uec.ac.jp

今回は、無料かつオフラインで使えるMaximaを解説してみます。いつの間にか僕の愛用ソフトになっていたんだ...

もともと難しい計算は wolfram alphaに投げていたのですが、応答が長かったり、変数を記憶させておくことができなかったりと不満点がありました。そこで、オフラインで使えるMaximaを主に使うようにしたという経緯があります。同じような不満を抱いている方にはこの記事は需要がありそう。

手計算って死ぬほど効率悪いですし、電通大生なら愛用する数式処理ソフトを一つ作っておくのは生き残り戦略的にアリです。計算ミス減らせるので、課題ほぼ満点狙えるのは強すぎる。GPAボクシングに突っ込め！！！

Maximaを導入する

Maximaの導入方法については、調べれば記事はいくらでも出てくるので割愛します。Maximaだけでも数式処理はできるのですが、wxMaximaというものを導入してそれから操作すると、かなり使いやすくなると思うので、これについても導入をお勧めします。(記事もwxMaximaから操作して動かしています)

Maxima,wxMaximaを共に導入し、wxMaximaを起動するとこんな感じの画面が出ます。

f:id:bokuroro:20211215233827p:plain — wxMaximaの初期画面

ここに数式を入力して処理していきます。試しに $(x+1)^7$ の展開でもしてみましょうか。手計算だとめんどそうですね。式を入力し、セミコロン+ENTERを押して式を評価してみると、次のようになりました。

f:id:bokuroro:20211218021731p:plain

%i+数字と書いてあるのが、コマンドで入力した部分、%o+数字と書いてあるのがそれと対応した出力ですね。%とすると前の出力を次の入力とすることができますから、それをexpand(展開)コマンドで処理すると次のようになるわけです。

変数や関数を扱うこともできます。次の例はaという関数を定義して、yで1回偏微分している様子を表しています。diffというコマンドが微分を表しているわけですね。

f:id:bokuroro:20211218021827p:plain

まぁ色々コマンドはあるので、調べてみるといいと思います。何もかもできるとは限りませんが、大抵のことならできると思いますよ。

初心者の方で、使い方を覚えるのがめんどくさい！！という方はwxMaximaのメニューバーに式の展開とか、方程式を解くといったコマンドが簡単に入力できるようになっているのでそれを使っていくのがおすすめです。

簡単な紹介が終わったので実際の課題をどう解くかといった様子をお見せしたいと思います。具体的な例があった方がわかりやすいと思うので。

課題は楽に解ける！！！

線形代数編

次の実対称行列$A$を適当な行列$P$を持って、$P^{-1}AP$のように対角化せよ。

$$A = \begin{bmatrix}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \\ \end{bmatrix}$$

また、$P$をグラムシュミットの正規直交化法により直交行列に直せ。

固有値とか固有ベクトルを計算するのとかってだるいですよね... 例えば行列Aの固有値が知りたかったら、$\lambda$についての方程式

$$|\lambda E - A| = 0$$

を解かなくてはいけないのでした。だるい！！　グラムシュミットの正規直行化法も知っている人ならわかると思いますが、めちゃくちゃだるいです。でもMaximaなら一瞬。

次のコマンドで行列Aを入力します。A:とすることで、計算結果をAという変数に代入することができます。

(%i)  A:matrix([2,1,1],[1,2,1],[1,1,2]);

その後、固有ベクトルを次のコマンドで計算します。

(%i) eigenvectors(%);

結果としては次の通りになりました。

f:id:bokuroro:20211217223145p:plain

わかりにくいですが、こんな感じです。

f:id:bokuroro:20211217223540j:plain

配列のような形で表示されていて、一番初めが固有値、次がそれぞれの固有値の重複度、その次がそれぞれに対応する固有ベクトルとなります。(ただし、行ベクトルとして出力されるので、後で転置を行う必要があることに注意)

固有値$\lambda = 4$(重複度1)の固有ベクトル

$$\begin{bmatrix}1 \\1 \\1 \\ \end{bmatrix}$$

固有値$\lambda = 1$(重複度2)の固有ベクトル

$$\begin{bmatrix}1 \\0 \\-1 \\ \end{bmatrix},\begin{bmatrix} 0 \\1 \\-1 \\\end{bmatrix}$$

となります。

よって、下のコマンドにより$P$を求めてみます。少し複雑ですがeigenvectorsの出力結果がリストになることに注意して読んでみるといいと思います。詳しい説明としては、一行目で固有ベクトルのみを抽出した後、二行目で行ベクトルをつなぎ合わせて行列Pを作成しています。列ベクトルに直したいので、行列Pの転置を行っているわけです。ややこしいなら行列$P$を手入力しても良いと思います。

(%i) [vecnum,vects]:eigenvectors(A);
(%i) P:matrix(vects[1][1],vects[2][1],vects[2][2]);
(%i) P:transpose(P);

結果としては次の通りです。$P^{-1}AP$も計算してみます。行列の積は .(ドット)、逆行列はinvertコマンドで表します。

f:id:bokuroro:20211217225528p:plain

よって、対角化行列が求められました。

グラムシュミットの正規直交化法

また、この行列Pをグラムシュミットの正規直直交法を用いて直交行列を求めてみます。グラムシュミットの正規直交化法は標準のMaximaの機能では計算できないので、eigenパッケージを導入します。

(%i) load(eigen);

このeigenパッケージは固有値や固有ベクトルの計算が得意なパッケージらしいです。(筆者調べ)

グラムシュミットの直交化法を行う(正規化はしてくれないので不親切)コマンドであるgramschmidtを実行します。

(%i) gvects:gramschmidt([vects[1][1],vects[2][1],vects[2][2]]);

これらをそれぞれ正規化(unitvectorコマンド)して、直交行列を求めると次のようになります。

(%i) transpose(matrix(unitvector(gvects[1]),unitvector(gvects[2]),unitvector(gvects[3])));

f:id:bokuroro:20211217234659p:plain

よって、直交行列$U$を得ることができました。一応、${}^t U A U$を求めてみましょう。

f:id:bokuroro:20211217234829p:plain

は？？？？？？？？？？？？？？

深夜に書いているのでつい怒りが出てしまいました。

数式処理ソフトといえども、式の展開や計算までは自動的にはやってくれないのです。人間の「行列の積を計算しろ！！」とか「因数分解しろ！！」*1とかいう指示に従って動いているだけということを常に意識しておくと良いです。勝手に人間に都合が良い形にしてくれるほど賢くはない。

この場合はratsimp(式の整理)とexpand(展開)をかければ良いです。

f:id:bokuroro:20211217235240p:plain

よかった、無事対角化できることがわかりましたね。

力学演習編

慣性モーメントを求めるのって重積分とかいきなり出てくるからかなりキツイですよね。次の簡単な例をもとに計算してみましょう。

密度$\rho$で一様で、半径$a$の球の中心を通る軸周りの慣性モーメント$I$を求めなさい。

力学演習取ってる方なら、慣性モーメントの定義は大丈夫ですよね？こんな感じになります。

$$I = \int_{D} \rho (x^2 + y^2) dx dy dz $$

ただし

$$D : 0 \leq x^2 + y^2 + z^2 \leq a^2$$

Maximaは重積分を計算できないので累次積分に直す必要があります。手計算だとどういったやり方をしていたでしょうか。

これを求めるには極座標変換をするのでした。

$$x = r \sin \theta \cos \phi, y = r \sin \theta \sin \phi, z = r \cos \theta$$

また、領域$D$は領域$E$に変換できて

$$E : 0 \leq r \leq a, 0 \leq \theta \leq \pi, 0 \leq \phi \leq 2\pi$$

ヤコビアンに注意して

$$dxdydz = r^2 \sin \theta dr d\theta d\phi$$

こんな感じです。準備が整ったのでMaximaで計算していきましょう。

被積分関数を定義し、極座標変換します。

(%i) f : rho* (x^2 + y^2);
(%i) g:f, x=r*sin(theta)*cos(phi),y=r*sin(theta)*sin(phi);
(%i) g:trigsimp(g);
(%i) J:r^2*sin(theta);

2行目の命令は、fの中のxとyをこの式に置き換えよ、という命令になっています。つまり、曲座標変換をしているわけですね。また、3行目の命令は三角関数の相互関係を使って式を整理するコマンドです。最後のはヤコビアンを定義しています。

f:id:bokuroro:20211218013357p:plain

よって、これより積分を実行していきます。定積分の計算にはintegrateコマンドを使います。

(%i) g*J;
(%i) integrate(%,r,0,a);
(%i) integrate(%,theta,0,%pi);
(%i) integrate(%,phi,0,2*%pi);

一行目は被積分関数とヤコビアンをかけて、二行目以降で積分を実行しています。二行目は$r$について、3行目,4行目は$\theta,\phi$についての積分になっています。($\pi$は%piと表すことに注意してください)

結果は次の通りです。

f:id:bokuroro:20211218013849p:plain

よって

$${{8\,\pi\,a^5\,\rho}\over{15}}$$

となりました。

これは簡単な例でしたが、他の重積分を求める問題に応用できるので、例えば工基礎の面積分の問題とかに応用してみてください！

工学基礎数学編

二年前期の工学基礎数学、計算量が多くて死にそうになってると思いますが、簡単に解いて楽しちゃいましょう！！

$$\boldsymbol{v} = (-x^3 y^2 - 4xyz^3) \boldsymbol{i} + x^2 y^3 \boldsymbol{j} + y z ^4 \boldsymbol{k}$$

の流れ(ベクトル場)について次の問に答えよ。

(1)流れが渦なし(rot = 0)かどうか、また非圧縮系(div = 0)であるかどうか調べよ。

(2)流れが渦なしなら速度ポテンシャル(スカラーポテンシャル)を、非圧縮系ならベクトルポテンシャルを一つ求めよ。

工基礎でよく見るベクトル解析の問題ですね。(僕が受けた問題から適当に引っ張ってきました。問題があれば消します。)rotとか二度と手計算で計算したくないですね。では計算していきましょうか。

ベクトル解析用のパッケージ、vectパッケージを呼び出します。

(%i) load(vect);

(1)

少しややこしいですが、次のように計算していきます。まず、ベクトル$\boldsymbol{v}$を入力します。

(%i) v:[-x^3*y^2-4*x*y*z^3,x^2*y^3,y*z^4];

まず、rotを計算します。Maximaは、rotのことをcurlと表しているので*2 、次のように計算できます。

(%i) curl(v);
(%i) express(%);
(%i) ev(%,diff)

実行結果としては以下の通りです。

f:id:bokuroro:20211218001152p:plain ややこしいですが、curlと書かれたコマンドを、expressコマンドで微分に変換し、evコマンドでdiff環境で式を評価し、微分を実行するという流れになっています。(ややこしいので深く理解しなくて良いと思います。こういうふうにするとrotが計算できると覚えちゃうのもアリ)