December 3, 2004

数式処理システム Maxima で材料力学

Maximaというのは、数式処理(Symbolic Computation)システムのアプリケーションである。中谷彰宏の数式処理には以前からずっと興味を持っていながらなかなか深いつきあいにならない。また、昔話になるが、学部のころmuMATH/muSIMPというソフトがCP/MというOSで動き、コンピューター上で数式処理ができることを知った。当時、恩師北川浩先生の友人の中村喜代次先生の研究室では(8bit から16bit への移行期で、おそらく8bit機との互換性を意識されていたのでは無いかと思うが)CP/M一色であり、先生どうし楽しく張り合っていたので、うちの研究室に配属になった時も当然8インチの大きなフロッピーのCP/Mのシステムディスクが転がっていた。ただし、大方はMS-DOSに移行していたように思う。当時、シンボリックな演算ができるmuMATH/muSIMP には感動したが、計算機のパワーを考えてもわかることだがあまり賢い計算はできなかった。M1のとき、Sony のワークステーション NEWS-831 を買ってもらったとき、ついでにREDUCE を入れてもらった。これは、そこそこ使えた。研究ではM2の時に、当時4年生の阪本正悟君に分子動力学シミュレーションの結果の結晶方位決定に着目原子の隣接原子の配置と参照配置の間での差の2乗ノルムのようなものを目的関数にし最小二乗法で最小化するというベタな方法を思い付いてよろこんでややこしい多変数関数の偏微分をやらせた。ニュートンラフソンのためのヘッシアンのFortran出力をそのままプログラムに貼ったら解読不可能のものすごい量になった。その他、杉山吉彦先生の弾性安定論の講義のレポートで平板翼のフラッター流速をHurwitzの方法で求める時に使わせてもらった(良く知らないけれど、この手の計算はMATLABでやるのが普通なのだと思う)。また、「数理科学」などにも数式処理のできる電卓の広告がいつも出ていたので欲しいなと思いながら「でも何に使うの?」ということで思い留まった(Brown Book Store on Thayer Street でもだいぶ迷った)。次は、MathCADで、これはモニターだったので、何かレポートを提出した。Maple は使ったことがなく、Mathematica もうちの大学でサイトライセンスをとっているから、導入の敷居はかなり低いことは知っているのだが、なかなか使う気にならなかった。そんな折、Maxima というのがフリーでLinux 上でも動く非常にしっかりしたソフトであることを知った。Richard Stallman 博士程ではないが日頃フリーソフトの重要性や意義を説いている自分としてはソースコードも公開されていることは重要である(とはいえ、Lispと名のつくのはemacs-lispぐらいしかみたことがないのでよくわかっていない)。事情を覚えていないがこれを知ったきっかけは、「Linux Japan」での松田七美男先生の連載だったのではと思う。それから、昨年の卒業生の小林佳広君と一緒にああだこうだ言いながら結構使えるなということになった。 それで、しっかり、ANC9(ANCollections Ver.9; 中谷彰宏が研究の現場で何かと便利なツール等集めて来たり、必要に応じて自前でrpmにbuildしたもの。まだredhat9版でとまっており fedora core に完全に移行できないでいる。余談になるが、研究室のマシンを全てFCに移行するふん切りがつかない理由はいくつかある。初期のころインストーラーがこけるマシンが多かったこと、ABAQUS, Tecplot などの商用ソフトを導入しこれらがFCで安定して動かせる保証、自信がなかったこと、などである。)に、maxima, maxima-exec-gcl, maxima-xmaxima, TeXmacs, TeXmacs-fonts などのパッケージを加えていたが、しばらく忘れていた。昨日、学生に材力の講義で黒板に式を板書しているときに、ふと頭をよぎって。。。講義を終るや否やすぐに試したかったが、いろいろとやっているうちに遅くなった。

材料力学関連の研究分野では、松本英治先生、片山忠一先生、菅野良弘先生のご研究の講演を聞いたようなおぼろげな記憶があるが最近の動向については全く不勉強である。 一方、材料力学関連の教育分野については、検索エンジン等で少し調べたところ、材力を数式処理で扱う話は、いろいろとあるらしい。当然だと思う。また、Mathematica を使うことを前提にした材力の本もあるとのこと。ただ、"Strength of Materials" や、"Structural Mechanics" と "Symbolic Computation" の組み合わせでは、興味あるサイトを見つけられなかったので(もっと真剣にやればぼろぼろ出てくるに違いないが)、まずは、思い付くままやってみることにした。

beam 例題として選んだのは、等分布荷重pを受ける長さl曲げ剛性EIの一端固定他端支持はりの問題である(図)。 支持点での支点との間の不静定内力(すなわち不静定反力)をXとおいて、最小仕事の原理(この場合は、カスチリアーノの定理で解いても同じ形)を使って、支持端での反力を求める。Maxima(実際には、xmaxima)を使った手順は、以下の通り。使っている機能は、未知変数を決定した後の代入や再定義などで、もっとスマートな方法があるかもしれないが、なんせマニュアルを見るにも概念がわかっていないので、要求する機能からの逆引きができない。

(C1) RATFAC:true;   部分因数分解をする。
(D1)                                 TRUE
(C2) define(F(x),-p*x+X);   せん断力分布を外力と不静定内力の項で表現する。
(D2)                            F(x) := X - p x
(C3) define(M(x),integrate(F(x),x));   曲げモーメント分布をせん断力の積分で定義する。
                                               2                   (問題から積分定数を決定)
                                            p x
(D3)                          M(x) := x X - ----
                                             2
(C4) AnsX:solve(integrate(M(x)*diff(M(x),X)/(E*I), x,0,l), X);   単位荷重の定理による
                                       3 l p                    不静定内力の評価。
(D4)                              [X = -----]
                                         8
(C5) subst(AnsX,F(x));define(F(x),%);   せん断力の再定義。
                                  3 l p
(D5)                              ----- - p x
                                    8
(C6) 
                                      3 l p
(D6)                          F(x) := ----- - p x
                                        8
(C7) subst(AnsX,M(x));define(M(x),%);   曲げモーメントの再定義。
                                             2
                                3 l p x   p x
(D7)                            ------- - ----
                                   8       2
(C8) 
                                                 2
                                    3 l p x   p x
(D8)                        M(x) := ------- - ----
                                       8       2
(C9) define(i(x),integrate(-M(x)/(E*I),x)+C[1]);   たわみ角曲線(積分定数を含む)。
                                    3          2
                                 p x    3 l p x
                                 ---- - --------
                                  6        16
(D9)                     I(x) := --------------- + C
                                       E I          1
(C10) define(v(x),integrate(i(x),x)+C[2]);   たわみ曲線(積分定数を含む)。
                                  4        3
                               p x    l p x
                               ---- - ------
                                24      16
(D10)                  v(x) := ------------- + C  x + C
                                    E I         1      2
(C11) AnsC:solve([subst([x=0], v(x)=0), subst([x=l], v(x)=0)],[C[1],C[2]]); 
                                                   境界条件を使って積分定数を決定。
                                     3
                                    l  p
(D11)                       [[C  = ------, C  = 0]]
                               1   48 E I   2
(C12) subst(AnsC,i(x));define(i(x),%);  たわみ角曲線の再定義。
                              3          2
                           p x    3 l p x
                           ---- - --------     3
                            6        16       l  p
(D12)                      --------------- + ------
                                 E I         48 E I
(C13) 
                                  3          2
                               p x    3 l p x
                               ---- - --------     3
                                6        16       l  p
(D13)                  I(x) := --------------- + ------
                                     E I         48 E I
(C14) subst(AnsC,v(x));define(v(x),%);   たわみ曲線の再定義。
                               4        3
                            p x    l p x
                            ---- - ------    3
                             24      16     l  p x
(D14)                       ------------- + ------
                                 E I        48 E I
(C15) 
                                   4        3
                                p x    l p x
                                ---- - ------    3
                                 24      16     l  p x
(D15)                   v(x) := ------------- + ------
                                     E I        48 E I
(C16) RAT(v(x));   少し綺麗な形に(factor(v(x));にして因数分解しても良いけど)。
                               3        2    3
                           (2 x  - 3 l x  + l ) p x
(D16)/R/                   ------------------------
                                    48 E I
(C17) quit();

紙と鉛筆とを全く使わずにここまでできるのは感動ものである。ちなみに、texmacs というのを小林佳広君が見つけてきたのだが、これを使うと、こんなふう( texmacsの出力のスクリーン例 )に数式が綺麗に出る。こちらもなかなかの感動ものだ。

また、以下のように簡単に、せん断力図(Shear Force Diagram, Shearing Force Diagram; SFD)、と曲げモーメント図(Bending Morment Diagram; BMD)やたわみ曲線(Deflection Curve)を書くことができる。

plot2d(subst([l=1,p=1], F(x)),[x,0,1]); 
plot2d(subst([l=1,p=1], M(x)),[x,0,1]); 
plot2d(-subst([l=1,p=1,E=1,I=1], v(x)),[x,0,1]);

あまりに感動して、落ちをつけられなくなったのだが最後に一言。 「材料力学応用」を受講している2年生と自分とで、この手の計算の能力がどちらが優れているかといえば、残念ながら彼らに軍配が上がる。昨日うちの研究室の学生とも雑談したけど、個人差はあるに違いないが、このような初等解析学や高速演算処理能力のピークは、入試の直後から大学院前期課程の入試までピークというのは間違いないように思う。はっきり言って平方根の展開など綺麗さっぱりわすれてしまっている(使っているCanna の辞書では「因数分解」が変換されない。ほとんど自分の中で死語になっている?)。このようなツールをありがたく思っているのは僕だけ?。

[2004/12/6追記] なお、はりの途中に集中荷重や集中モーメントを受ける場合は、区分定義関数(piecewise defined function)が扱えなければならない。IF 〜 then 〜 else 〜 で関数の定義はできるものの、その微分、積分となるとうまくいかないようだ。これぐらいはできて欲しいが、これについては、assume で、変数の定義域を連続部分内で規定してやると良いようだが、あまり便利でないように思われる。ヘビサイドのステップ関数なども余計にややこしくなりそうで避けたい。もっとも、手計算でも同じように場合わけしているので繁雑さは似たり寄ったりではある。

[2005/01/19追記] 変わり映えしないけど、Maxima で材料力学(第2弾)もやってみました。

投稿者 Akihiro Nakatani : December 3, 2004 3:43 PM