模造２（偏微分方程式）
このページはマス旅の一部です。
今回は「偏微分方程式」でシミュレーションしよう。
まず、イメージを作ろう。
＜一般的なイメージ＞
常微分方程式（ODE)は1変数の微分方程式です。
geogebraでもそうですが、プログラミング言語には、数値的に解決する関数が用意されています。
解析的に解くと大変でもなんとかなりましたね。
しかし、偏微分方程式（PDE)はそう簡単ではありません。
2変数になると同じ方式が使えないのです。∂^2U/∂x^2をUxxなどとかくことにしよう。
2階偏微分方程式に限ると3種類に分類できます。
・電磁気など場の関数を扱うのが「ポアソン方程式（Uxx+Uyy=f(x,y)）」です。
その特殊型でf=0を特にラプラス方程式と言いますね。
式の形から楕円型というラベルを貼るときもあるね。
・時間の関数としての波動（音や光など）を扱うのが波動方程式Utt=cUxxです。
式の形から双曲型と呼ばれます。
・熱や物質の拡散を扱うのが拡散方程式Ut=cUxx。
これは放物型と呼ばれます。
常微分方程式では初期値で曲線関数y=f(x),y=f(t)などが決定しました。
偏微分方程式では境界条件で曲面関数、曲線関数z=f(x,y),U=f(x,t)などが決定するでしょう。
＜解析解のイメージ＞
たとえば次のようなステップを踏みます。
2階2変数の偏微分方程式を、
U(x,y)=g(x)h(y)の置き換えと微積（変数分離と部分積分など）によって変数分離する。
2階常微分方程式を2本できる。
別々に一般解の数式を求めて、境界条件から特殊形にする。
U(x,y)=Σg(x)h(y)[n for 1 to ∞]の形で合体する。
フーリエ級数を使って、f(t)に置き換える。
境界条件をfにも使うことで、最終的な曲面Uの方程式が関数の総和
でえられます。くわしくはこちら（微分方程式の境界値問題を解く）
常微分方程式を解くモジュールがあっても、
その前の分解作業は数式の特徴に合わせた処理が必要になります。
また、最後の合体したあとのフーリエ級数の活用でも
三角関数と指数関数に絡む積分と総和処理が絡むので、モジュールに投げることはできません。
＜格子点の差分法のイメージ＞
数値計算では、ラプラス方程式を例にとります。
２階微分が０なので、局所的にまっすぐな曲面です。
曲面の高さを数値に置き換えてイメージしよう。
１から100までの数を10個で改行してかくと、10×10の正方形の
数表ができますよね。その数表で十字型に数を5個選びましょう。
１つの数の上下左右に数をくっつけた十字型です。
十字の数は必ずたて、よこ、ななめに等差数列になるので、
中央が平均になりますね。等差数列というのは直線の証拠です。
これがラプラス方程式でも同じです。
Uxx+Uyy=f(x,y)の解の曲面がU(x,y)だとします。
境界条件がU(x,0)=0,U(x,1)=x,U(0,y)=0,U(1,y)=yのとき、
x,yが0以上1以下の正方形を領域Dと呼ぼう。
本当はDを微細な切片に切り刻むのですが、
たとえば、4等分線をたてよこ引くと、線分の両端を含むことで、
格子点は5×5＝25個できるね。
その高さUを離散的な数値として、行列にすることができる。
境界線を含む点は境界条件から自明だが、それ以外の
内部点3×3＝９個は自明ではない。この9点を行列にできる。
A={{u11,u12,u13},{u21,u22,u23},{u31,u32,u33}}
ここで、中央が平均だという条件から、
9個の点の高さuijが他の4つの変数や数値の平均の式に直せる。
だから、9元連立方程式になるね。
これを行列にするとAx=bの形でAが9×9の正方行列ができる。
これを解くためには[A|b]の行列をガウス消去法でカンタンに直せば、順次的に９つの元が出せる。
ということは境界以外の内部の曲面の高さが決まるということだね。
課題１：熱の伝導（拡散方程式・放物型）
「f(x)=if(x≦L/2, x, L-x)とする。
（境界条件）x=0, x=Lでは温度が常に０に保たれている。
u(0,t) = u(L,t) = 0
（初期条件）ｘ＝L/２で最高になる三角形の分布。
f(x)=u(x,0)
f(x,t)の関数を決定して時間ととも温度分布がどう変わるかを視覚化しよう。」
＜解析解＞
そこで、上のイメージ通りの手順でやると、
u(x,t) =
nが偶数ならCn=0,
n=1,5,9,....,ならCn , n=3,7,....なら、Cn
＜コード化＞
n=奇数の項をCnの符号を変えて出し、それをSumすることでu(x,t)のグラフになる。
実際はtをスライダーにして、uが縦軸、xが横軸のグラフにするといいね。
L=π^2などとおけば、数式は簡単になるでしょう。
#PC環境によってはフリーズしたようになるかもしれません。
関数列の個数を10個までにしてますが、環境次第では5とかに相当下げてください。
プレイボタンはストップにしてあります。
// 1. 基本パラメータとスライダーの設定
L = pi^2
c = 1
t = Slider(0, 5, 0.1) // 設定 ＞ 上級 ＞ アニメーション ＞ 増加 に変更
// 2. 係数リスト Cs の作成 (k=1,3,5,...の奇数項のみ値を持ち、符号が交互に変わる)
Cs = Sequence(If(Mod(k, 2) == 1, (-1)^((k - 1) / 2) * 4 * L / (k^2 * pi^2), 0), k, 1, 10)
// 3. 各項の関数をつくり、Sumで総和をとる
u_t = Sequence(Element(Cs, k) * exp(-(k * pi / L)^2 * c * t) * sin(k * pi / L * x), k, 1, 10)
u(x) = if(0<=x<=L,Sum(u_t),0)
中央最大化の角張った温度分布が
時間ととものまろやかになるのが曲線の形から「見たまんま」わかるね。
拡散方程式
課題2：弦の振動（波動方程式・双曲型）
「u=u(x,y)について utt = c^2 uxx
初期条件 初期変位u(x,0)=f(x) 、初期速度ut(x,0)=v(x)。
初期変位f(x)=if(x≦L/2, x, L-x)とする。
境界条件u(0,t)=u(L,t)=0 から、f(0)=f(L)=0となる.
弦の振動を視覚化しましょう。」
＜解析解＞
u(x,t)=Σ Cn cos(ncπ/L t) sin (nπ/L x)
x=0, x=Lは固定されていて
変位はｘ＝L/２で最高になる三角形の分布f(x)=u(x,0)
Cnは熱伝導と同じ式になる。
熱伝導のexp(-(nπ/L)^2 ct)部分が
波動ではcos( (nπ/L) ct)に置換されたともいえるね。
スピード
＜コード化＞
// 1. パラメータとスライダー
L = pi^2
c = 1
t = Slider(0, 5, 0.05) // アニメーション：増加
// 2. 符号リストと係数列の作成（定数部分は外に出す）
sg = {1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0}
cr = Sequence(sg(k) / k^2, k, 1, 9)
// 3. 各項の段階で If を組み込み、Sum で一気に合成
// （※ 係数 4L / pi^2 は最後にまとめて掛けるとさらに計算が軽くなります）
fn = Sequence(If(0 <= x <= L, (4 * L / pi^2) * cr(n) * cos(n * c * (pi / L) * t) * sin(n * (pi / L) * x), 0), n, 1, 9)
u(x) = Sum(fn)
中央最大化の角の引っ張りが、波として両端に伝わるようすがわかるね。
両端で跳ね返りながら減衰せずに繰り返すのがわかるね。
波動方程式
課題３：ゴム膜の収縮（ラプラス方程式・楕円型）
「z=U(x,y)についてUxx+Uyy=0
境界条件がU(x,0)=0,U(x,1)=x,U(0,y)=0,U(1,y)=y
となるピンと張ったゴム膜Uをかこう。」
＜格子点と連立方程式＞
D=[0,1]×[0,1]の辺の4等分点の3×3=9個の交点のzを3行3列の要素u_ijの行列M={u_ij}にする。
格子点u_ijどうしの十字平均の結果から、9要素a_klの連立方程式を作る。
連立方程式を行列とベクトルでAx＝bとかける。これを[A|b]についてガウス消去法を使い順次解けるようにする。[A|b]⇒[B|b']のとき、Bは主対角成分が1の右三角行列になる。ベクトル b'を利用して、最下行から消去法でｘの上位成分を順次求めることができ、ベクトルxとなる。
このｘが9要素がMの値をフラット化したものとなる。
答えは
{0.0625, 0.125, 0.1875, 0.125, 0.25, 0.375, 0.1875, 0.375, 0.5625}
これをN+1等分点のN*N個の交点について、同じロジックで連立方程式を作り、ガウス消去法などで方程式を解けば、領域内の各点(x,y)の値zが離散的にわかる。これらの点を通る面がかければよいね。
たとえば、x=cで点(c,y,z)を折れ線で結ぶ。このcを動かす。同様にy=dでxを動かして折れ線で結ぶ。
この独立作業によって、空間内をメッシュすることができる。
あるいは、4点指定のポリゴンリストを作ることで平面の断片をつなぐことできるでしょう。
連立方程式を解くにはSolveか行列のReducedRowEchelonFormコマンドが使えるでしょう。
しかし、多変数の連立方程式をそもそも組み立てたり、解いたりするのはgeogebraでは荷が重い。
そこで、「平均反復法」という裏技を使おう。
十字の「5点平均」を順次求めると、だんだんと平坦になっていく。
つまりラブラス方程式に迫れるということだね。
＜コード化＞
# 立体図形の作成
#解析解の面の表示
showExact = true
U_{exact}(x, y) = If(0 <= x <= 1 && 0 <= y <= 1, x * y, undefined)
ExactSurface = Surface(x, y, U_{exact}(x, y), x, 0, 1, y, 0, 1) #水色で濃さ50％#表示条件showExact
//内部点 3x3 の高さリスト U (初期値)
U = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}
// 境界条件を含む 5x5 全体の高さリスト Z
// Index 1..5 (0, 0.25, 0.5, 0.75, 1.0)
Z = {{0, 0, 0, 0, 0}, {0, Element(Element(U,1),1), Element(Element(U,1),2), Element(Element(U,1),3), 0.25}, {0, Element(Element(U,2),1), Element(Element(U,2),2), Element(Element(U,2),3), 0.50}, {0, Element(Element(U,3),1), Element(Element(U,3),2), Element(Element(U,3),3), 0.75}, {0, 0.25, 0.50, 0.75, 1.00}}
// 3D格子点およびメッシュ（MeshX,MexshYを表示にする。色は黒で太さ３）
GridPoints = Sequence(Sequence(Point({(i-1)/4, (j-1)/4, Element(Element(Z, j), i)}), i, 1, 5), j, 1, 5)
MeshX = Sequence(PolyLine(Element(GridPoints, j)), j, 1, 5)
GridTransposed = Sequence(Sequence(Element(Element(GridPoints, j), i), j, 1, 5), i, 1, 5)
MeshY = Sequence(PolyLine(Element(GridTransposed, i)), i, 1, 5)
#2Dのビューを表示にする。#「平均ボタン」と「リセットボタン」を貼り付ける。
#「平均ボタン」を押したときgeogebraScriptは以下を貼り付ける。見出しは「5点平均」
u11 = (0 + Element(Element(U,1),2) + 0 + Element(Element(U,2),1)) / 4
u12 = (Element(Element(U,1),1) + Element(Element(U,1),3) + 0.25/4 + Element(Element(U,2),2)) / 4
u13 = (Element(Element(U,1),2) + 0 + 0.50/4 + Element(Element(U,2),3)) / 4
u21 = (0 + Element(Element(U,2),2) + Element(Element(U,1),1) + Element(Element(U,3),1)) / 4
u22 = (Element(Element(U,2),1) + Element(Element(U,2),3) + Element(Element(U,1),2) + Element(Element(U,3),2)) / 4
u23 = (Element(Element(U,2),2) + 0.50 + Element(Element(U,1),3) + Element(Element(U,3),3)) / 4
u31 = (0 + Element(Element(U,3),2) + Element(Element(U,2),1) + 0.25/4) / 4
u32 = (Element(Element(U,3),1) + Element(Element(U,3),3) + Element(Element(U,2),2) + 0.50/4) / 4
u33 = (Element(Element(U,3),2) + 0.75 + Element(Element(U,2),3) + 0.75/4) / 4
SetValue(U, {{u11, u12, u13}, {u21, u22, u23}, {u31, u32, u33}})
#「リセットボタン」を押したときgeogebraScriptは以下を貼り付ける。
SetValue(U, {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}})
3軸とパースペクティブの自動計算の関係で、全貌が表示されるサイズはとても小さいです。
それでも、回転して眺めると、「5点平均」を押すたびにｘｙ平面からUのメッシュが解析解に
徐々に接近するのがわかりますね。「showExact」のオン/オフ切り替えボダンを使うと、解析解の表示/非表示が切り替えられます。