波動方程式シミュレーションをエクセルの循環参照と3Dグラフで。

以前の日記で、1次元空間(2次元時空)での波動方程式の数値解をExcelで表現しましたが
今度は次数を1つ上げて、2次元空間(3次元時空)での波動方程式の数値解をExcelで表現してみようと工夫しました＾＾


 
デモ画像はこちらになります。

まあマックスウェル型(古典型)、といいますか非量子・非シュレディンガータイプの波動方程式といいますか
いわゆるfutooの波が従う波動方程式なので、縦軸を変位とかと捉えれば3次元空間、4次元時空と考えられなくもない・・・


10/6のときと同様、循環参照を使って逐次計算(時間領域)を行っています。
前回は横：位置、縦：時間だったのですが、
今回は横も縦も位置で、シートに時間を対応させています。
 
循環参照を使うので、
・初期0と1　(1)
と、
・初期を参照する2と3　(2)、
そしてさらに
・初期(1)に頼らず(2)を参照する4と5　(3)
の時間領域を使います。

初期条件(1)を離れると、2、3(2)→4、5(3)→2、3(2)→4、5(3)・・・といった風に紡がれます。
 
これを今回は縦ではなくシートごとに行うことで次数を上げました。


＝＝＝＝＝＝＝＝＝

境界条件

また、次数が上がると境界条件の数もぐんと増えます。
時間に対しても位置に対しても2階の微分方程式なので、
前回の1次元空間＋1次元時間の場合は2つの初期条件と2つの境界条件でした。

これが今回次数を上げると、時間の初期条件が2つは同様ですが
空間の境界条件は4つになります。
 
 
（厳密には4面と4本と8点、でしょうか？）


今回はすべての境界においてロビン境界条件を一様に設定しましたが
 x、yの範囲や座標軸の違いなどでディリクレ条件にするかノイマン条件にするか選べるようにしてもいいと思います。(これが混合条件です。ロビン境界条件とは異なります)



ロビン条件を設定したとは言いましたが、実質的にはディリクレとして使ってます。
これを空間的に端っこまで表示しないことで、擬似的なロビン条件を満たします。
(インピーダンスマッチングがとれているので反射がない状態)



 
シミュレーションとして動かす際は、適当な時間(シート)・空間(セル)を指定して、そのセルに0でない数値を入れてください。
もし既に何かしらの式が入っている場合は、色を変えるなりして「ルールを変更した目印」をつけておいて、あとから上下左右の隣からコピペして修復しやすいようにしてください。



＝＝＝＝＝＝＝＝＝

循環参照

今回は0シートのy=0(10行目)のx=10とx=40(それぞれO列とAS列)に6という数値を入れてみました。
こうすることで、憧れの二重スリット実験のような波紋ができることになります。
動かす際は「循環参照をオンにして」スイッチ(sw)を入れて初期条件にリセットしたのち
空白セルでdelボタンを押しっぱなし(再計算)にしてください。
もしノロマだと思ったら、循環参照のオプション変更で、「最大反復回数」を調整してみてください。


 
主な操作はシート「5」でだいたい行えるようにしてあります。
 クーラン数(cdt/dx)^2は特に弄ることはないと思います。クーラン数が大きすぎると発散します。


＝＝＝＝＝＝＝＝＝

DL版のいじりかた(カッペーではなくコピペー)

DL版Excelファイルを用意しておきました。
DL版は容量の関係で、縦と横の領域幅が狭いです。
増やしたい場合は、「0シートをクリックしてからシフトを押しながら5シートをクリックしてすべてのシートの選択」をし、一律に並行作業をさせます。 (シート一括操作：※通常時は解除してください)
 選択したすべてのシートが白くなります
させる内容はまず横を広げるほうからやりますと、
 「0から5までのシートが選択されていることを確認した上で」
 
「右端(境界条件)の列全体を広げたい先までカッペーではなくコピペー」
     
それから、新しい端っこの1歩手前まで一般領域をコピペして広げればOKです。

縦を広げる方法も同様です。
 

   
下端(境界条件)の行全体を広げたい先まで カッペではなくコピペです。


なぜカットではなくコピーしてからペーストなのかといいますと
ポインタ・・・という言葉を今なお使っていいのかわかりませんが
コピペだと中の人が変わるのに対して、カッペだと中の人が固定なのです。
つまり、カットアンドペーストしてから色々処理を行ってそのセルが死んでしまうと、
参照ができなくなって後々の作業妨害になってしまうのです。非常に作業効率が悪くなります。




＝＝＝＝＝＝＝

ナイキスト

今回のシミュレーションでは、なんかミョーに波がざわついてしまいました。
たぶんですが、これは領域が有限であることと、離散化されていることと、
波動方程式という周期化しやすいシチュエーションが相互に影響しあっているんだと思います。
つまり、時間・空間どちらにおいても、実空間(フーリエ変換前)・(周)波数空間(フーリエ変換後)ともに
周期化・離散化されているので
 適当に間引きしないと本質的にサンプリング定理に抵触してしまい、意味のないノイズが出てしまうんだと思います。たぶん・・・


かといって球(円？)面波でデータを間引きするというのも正直、無理数から整数を取り出すみたいでどうやったらいいかわかりませんし・・・まあとりあえずそっとしておこうということにしました。＾＾；
領域が2次元になりましたからねえ・・・



あ、それと、波が2次元に広がったことで、振幅が距離の1乗に反比例して小さくなっていることが(なんとなく)わかると思います。1次元空間では空間的な減衰はしませんでしたからね。




＝＝＝＝＝＝＝＝

差分法の陽解法

どのような演算をしているのかといいますと。
波の振幅をuという関数とします。
x方向の番号をi、y方向の番号をj、時間方向の番号をkとしますと

一般領域では、波動方程式
 

の偏微分を差分に置き換えて

Δx=Δyとしてから
i、j、k番目のuの式に変型しただけなんです。
 

位置の境界条件も案外さらっと行きまして
問題は2つ目の初期条件ですね。

どうも、テイラー展開を使って
 
こういう風にしてるみたいなんですよ。

にほんブログ村