目次

python, vpython と sf による偏微分方程式の数値解法

はじめに

任意の偏微分方程式を数値的に解く solveLaplacian(..)/solvePDE(.) 関数/itrSlvPDE(..) イタレータ関数を python 言語で作りました。その使い方、応用例を以下に示しますので御評価願います。

solveLaplacian(..)/solvePDE(..)/itrSlvPDE(..) 関数の概要

掲題の三つの関数は全て境界条件行列を与えて、偏微分方程式に対応する漸化式を繰り返し適用して数値解を求める python 関数です。

      ┌──────┐    +--------+
      │境界条件行列│    | 漸化式 |
      └───┬──┘    +---------
              ↓            ↓      
       +--------------------------+ 
       |    soveLaplacian(..)     | 
       | or solvePDE(.. )         | 
       | or itrSlvPDE(..)         | 
       +--------------------------+ 
                    ↓      
           dctRslt:数値解行列
                    ↓      
        make3dRowCol(dctRslt):数値解の可視化

python モジュール sf.py にまとめてあります。計算結果を 3D グラフとして表示する関数 make3dRowCol(..) も用意してあります。

solveLaplacian(..) の使い方の概要

ブロック実行

本稿では、下に示すような //@@ と //@@@ で囲んだコードの記述を何度も使います。これは //@@ と //@@@ で囲んだ部分の文字列を \#####.### と名付けたファイルにし、//@@@ の次に続く // で始まる行の文字列をコンソール・コマンド文字列として実行するという意味です。下の例では //@@ から //@@@ の範囲の文字列を \#####.### ファイルにし、そのファイルをカレント・ディレクトリの a.cpp にコピーし、a.cpp を cl.exe でコンパイルする処理を意味します。

//@@
    ・
    ・
//@@@
//copy \#####.### a.cpp
//cl a.cpp -I.\

私はこのような処理をエディタで自動的に行うようにしています。wz ver4.0 エディタでのマクロで実装しています。ctrl O + E で、この一連の処理を行わせています。GPL 条件で公開していますので、wz エディタを使える方はご利用ください

python のときは、実行オプションが付かないことが大部分なので「 //@@ ... //@@@ 文字列を copy \#####.### temp.py と copy させ、カレント・ディレクトリで python.exe temp.py を実行させる」ことまでを ctrl O + P に割り振ったエディタ・マクロで行わせています。sf でのブロック式を計算させるときは 「temp.se を作り sf.exe @@temp.se を実行させる」ことまでを ctrl O + S に割り振ったエディタ・マクロで行わせています。//copy... や //python temp.py などを書かずに //@@...//@@@ だけで済ませています。

以下では //@@ と //@@@ で囲んだプログラム・コードが何度も出てきますが、このような意味であることを御承知ください。

もっとも単純な solveLaplacian(..) の具体的な使い方を見てみましょう。solveLaplacian(..) 側に偏微分方程式を解く漸化式が内蔵されているので、下のように境界条件行列を与えるだけでラプラス方程式が解けてしまいます。偏微分方程式を解くコードとしては限界に近いコンパクトさだと自負しています。如何でしょうか。

//@@
import sf
N = 30
#N =  8
dctBoundary={}
for i in range(N):
    for j in range(N):
        dctBoundary[i,j]= True

# set center boundary value                             #     ----  
dctBoundary[N/2-1,N/2-1]=2; dctBoundary[N/2-1,N/2]=2    #   |  22  |
dctBoundary[N/2,N/2-1]=2;     dctBoundary[N/2,N/2]=2    #   |  22  |
                                                        #     ----  
# solve Partial Differential Eqation
dctRslt = sf.solveLaplacian(dctBoundary)
sf.make3dRowCol(dctRslt)                                # visualize the result
//@@@

上の小さな pthon コードを記述するだけで、境界条件に従うラプラス方程式の数値解を求め、 下のような、3D 表示まで行います。この 3D 表示はマウス操作により回転や拡大もできます。

このコードの dctBoundary 変数が境界条件行列です。4 行目から 10 行目までで、下のような行列値を設定しています。下の True は python での bool 値です。

True True True True True True True True  
                        
True True True True True True True True  
                                         
True True True True True True True True  
                                         
True True True   2    2  True True True  
                                         
True True True   2    2  True True True  
                                         
True True True True True True True True  
                                         
True True True True True True True True  
                        
True True True True True True True True  
    
    8x8 境界条件行列 dctBoundary

上の行列値で True は偏微分方程式を解くために漸化式を繰り返して適用する領域であることを意味します。真ん中の数値 2 は初期条件/境界値です。

ラプラス方程式での漸化式を適用できない端面での値は線形外挿補間で求めています。線形外挿補間の位置・方向は、偏微分方程式の求解領域が直方体であることを利用して、境界条件行列引数から solveLaplacian(..) 側で推測させています。これにより、 sf.py による偏微分方程式の数値解コードを限界近くまで簡素化できるようになりました。

solvePDE(..) の使い方の概要

ラプラス以外の一般の偏微分方程式を数値的に解くときは、偏微分方程式に対応した漸化式をユーザーが用意します。あえて上のラプラス方程式を solvePDE(..) を使って漸化式を明示したコードにすると下のようになります。fncPDE(..) 関数がラプラス方程式に対応した漸化式です。計算結果は solveLaplacian(..) のときと全く同じです。

//@@
import sf
N = 30
#N =  8
dctBoundary={}
for i in range(N):
    for j in range(N):
        dctBoundary[i,j]= True

# set center boundary value                             #     ----  
dctBoundary[N/2-1,N/2-1]=2; dctBoundary[N/2-1,N/2]=2    #   |  22  |
dctBoundary[N/2,N/2-1]=2;     dctBoundary[N/2,N/2]=2    #   |  22  |
                                                        #     ----  
def fncPDE(dctRslt, index):
    j,k=index
    dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                          + dctRslt[j+1,k]
                          + dctRslt[j,k-1]
                          + dctRslt[j,k+1])

# solve Partial Differential Eqation
dctRslt = sf.solvePDE(dctBoundary, fncPDE)
sf.make3dRowCol(dctRslt)                                # visualize the result
//@@@

波動偏微分方程式を解きたいときは fncPDE(..) 以下の漸化式を波動方程式に対応したものに変更します。漸化式さえ用意できれば、任意の偏微分方程式を解けます。

Python 言語での記述だけを考えるならば、次の itrSlvPDE(..) 関数を使ったコードの方が少しばかりきれいです。でも C 言語での高速化を後で行うことを考えているならば、solvePDE(..) でコードを書いたおいた方が、python から C 言語への変換の段階で楽になります。

itrSlvPDE(..) の使い方の概要

solvPDE(..) の実装の仕方は python での iterator 関数として実装したほうが、python 的です。記述も少しですが楽になります。

solvePDE(..) は境界条件に従って漸化式を繰り返し実行させています。漸化式関数は call back 関数として solvePDE(..) に引き渡しています。C 言語などで良く行われる実装のしかたです。

でも python に慣れた方ならばイタレータ:itrSlvPDE(..) を使って、繰り返し処理を下のように記述するほうを好むと思います。ここでも最初のラプラス方程式を例題に使います。計算結果も solveLaplacian(..) のときと全く同じです。

//@@
import sf
N = 30
#N =  8
dctBoundary={}
for i in range(N):
    for j in range(N):
        dctBoundary[i,j]= True

# set center boundary value                             #     ----  
dctBoundary[N/2-1,N/2-1]=2; dctBoundary[N/2-1,N/2]=2    #   |  22  |
dctBoundary[N/2,N/2-1]=2;     dctBoundary[N/2,N/2]=2    #   |  22  |
                                                        #     ----  
# solve Partial Differential Eqation
for dctRslt, index in sf.itrSlvPDE(dctBoundary):        
    j,k=index
    dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                          + dctRslt[j+1,k]
                          + dctRslt[j,k-1]
                          + dctRslt[j,k+1])

sf.make3dRowCol(dctRslt)                                # visualize the result
//@@@

solvePDE(..) と itrSlvPDE(..) の機能は全く同じです。繰り返し処理の記述の仕方が変わるだけです。ユーザーの目的・好みに応じて使い分けてください。

make3dRowCol(..) vpython:計算結果の表示

通常、偏微分方程式の数値解は 100 以上の要素からなる数値の塊りにすぎません。この数値の塊りだけでは、偏微分方程式の数値解を人間が理解することはできません。視覚化が必要です。

そのため sf.py モジュールにある make3dRowCol(dctRslt) 関数を使って、偏微分方程式のメッシュ行列数値解を 3D 表示させます。dctRslt 引数は、solvePDE(..) 関数が計算して返した偏微分方程式のメッシュ上の値を表す数値行列です。この 3D 表示のために vpython を使います。特に二変数の偏微分方程式までならば、ユーザーは vpython のコードを記述する必要はありません。make3dRowCol(..) 関数が dctRslt 行列の形から表示領域を自動的に定めて 3D 表示のために必要な vpython 処理を自動的に行います。

三変数以上の偏微分方程式の数値テンソル解については共通する視覚化ルーチンを作れません。単純に扱うと四次元以上の空間表示になってしまうためです。個別に視覚化コードを書かねばならない場合も多いでしょう。自分で vpython コードを書くことも多くなると思います。でも 3D 表示機能の塊である vpython パッケージを使えば、それも簡単です。ここでは vpython の

方法を探す必要があります。この原稿での応用例は二変数の偏微分方程式に限定します。でも soflvPDE(..) 関数自体は任意数の変数で記述される偏微分方程式を扱えます。div E のようなベクトル値の偏微分方程式さえ扱えます。

なお三階以上のメッシュ・テンソルであっても、特定の軸にそったベクトル数値ならば簡単にグラフ表示できます。特定の面にそった行列数値ならば 3D 表示できます。二次元 3D 表示にしても、特定位置での断面をグラフ表示させたいことも多くあります。このような処理をさせるとき、次に述べる sf が便利です。

python の必要性

汎用的な偏微分方程式の solver プログラムとして実装するためには python 言語を使わざるをえませんでした。solveLaplacian(..)/solvePDE(..)/itrSlvPDE(..) これらの関数が扱えるのはは二変数の偏微分方程式に限りません。N 変数の偏微分方程式を扱えます。ベクトル値の偏微分方程式さえ扱えます。このような汎用のプログラムの記述はは C などでは無理です。現実に存在していません。

計算結果の 3D グラフィックス表示のためにも python を使わざるを得ませんでした。make3dRowCol(..) は 150 行弱の python 関数です。必要ならばユーザーが自分で手を入れられる規模です。C などで同等の機能を実装しようとしたら一桁以上の大きなプログラムになるでしょう。ユーザーにも表示ライブラリの使い方を理解してもらうだけで、偏微分方程式の解法で使う以上の勉強時間が必要になってしまうでしょう。

上の solveLapracian(..) などのコードを見てもらえば分るように、C 言語などの素養があれば、python を使うのは簡単です。偏微分方程式を解くような決まりきった処理に限定すれば、後にも提示する幾つかの python コードの一部を修正するだけで多くの場合に対処できます。後で説明する「偏微分方程式を解くめの python 入門」を読む程度で、ユーザーが解こうとする偏微分方程式の漸化式を扱うことができます。

Python には scipy という数値計算に特化した無償のパッケージが存在します。行列計算、微分, 積分、フーリエ変換などが簡単に行えます。この機会に python を学んでおいても損をしないはずです。

C++ などを使えば計算速度が二桁近く速くなるのは分っています。しかし python でも 40x40 程度のメッシュならば一秒程度でグラフ表示でできるので十分に実用的です。python での記述のしやすさを考えると python のほうが勝ります。100 x 100 x 100 メッシュのような規模となると python 処理には向きませんが、小規模の試行錯誤の段階では python を強く勧めます。

sf:計算結果の検討

偏微分方程式を解いた結果としてできあがる膨大な数値の塊は 3D 表示だけでは把握しきれないことが多くあります。そのときに sf を使います。

sf は私の開発した行列・ベクトル計算を行わせるソフトです。エディタで書いている数式のままに行列ベクトル計算させます。これを使って、膨大なメッシュ行列/テンソル数値データから特定の直線/平面にそったデータを抜き出すときに、行列ベクトルデータの抜き出しや加減乗除算で、また それらの数値のグラフ表示で活躍します。

sf.py による偏微分方程式 数値解析の勧め

偏微分方程式の数値解を求めることで偏微分方程式が表現している自然法則を深く理解できるようになります。

Maxwell 方程式をはじめ多くの自然法則が偏微分方程式の形で記述されます。でも偏微分方程式を眺めているだけでは、その方程式が記述している内容を殆ど読みとれません。解析的に扱おうにも偏微分方程式を扱う統一的な方法、たとえば常微分方程式における Laplace/Fourier 変換のようなものはできていません。

でも現在のコンピュータ・パワーを利用すれば偏微分方程式であっても数値解析できます。個別具体例について偏微分方程式の数値解を求められます。実際に自分で具体例について、偏微分方程式の数値解をとき、それを様々の角度から視覚化することで、偏微分方程式が表現している自然法則が分ってきます。

でも偏微分方程式の数値解のプログラムの作成・デバッグに何時間何日もかけていたのでは本末転倒です。偏微分方程式で定まる物理・工学法則の理解を目的とするならば、そのためのプログラム・デバッグ作業は最小限のものでなければなりません。

solveLaplacian(..)/solvePDE(..)/itrSlvPDE(..) を使うことで漸化式を使った偏微分方程式の数値解プログラムの作成の手間を限界近くまで少なくできます。これらの関数では自由端の線形外挿補間を自動で行います。計算結果の可視化は、これら関数の戻り値行列値を sf.make3dRowCol(..) 関数に引き渡すだけです。ユーザーが用意するのは偏微分方程式の初期条件行列と漸化式の処理プログラムだけです。

Python 言語などの知識が必要な点で最初のハードルは高いかもしれません。でも その努力の見返りは十二分にあります。ぜひとも soveLaplacian(..)/solvePDE(..)/itrSlvPDE(..) に挑戦してみてください。 以下 python 言語の解説から始めていきます。

偏微分方程式を解くための python 入門

ここでの偏微分方程式を扱うために必要な最小限の python 言語についての解説を行います。C または Java 言語の素養があることを前提に話を進めます。

なお、以下での python 言語の説明は solveLaplacian(..)/solvePDE(..)/itrSlvPDE(..) を扱うコード例での一部を書き換えて、読者が別のコードを実行させられるようになることに目的として説明しています。Python 言語自体を理解することは、このような短い説明だけでは無理があります。また読者の言語素養によっては、下の python 言語の説明では理解してもらえないことも多いと思います。Python の知識がないときは、実際に python をインストールして自分の手で試しながら理解していく作業が必須だと思います。私の説明では不十分なときは自分でテスト・コードを様々に変形して試してみることが、結局は理解の早道だと思います。

python のインストールについて

Python2.5 は 2006年 9 月に出ています。でも今でも敢えて python2.4 の Enthought パッケージをインストールすることを進めます。多くの python モジュールが互いに矛盾することなくインストールされるからです。Python2.5 にするのは Enthought パッケージが 2.5 対応になってからで十分だと主張します。

Windows 環境の方は ここ の一番下にある enthon-python2.4-1.0.0.exe をダウンロードして実行してください。

python の文法

Python は C や Java 言語に似た手続き型の script 言語です。Script 言語の常として型宣言はありませんが、bool 値:True/False, 整数:int, 浮動小数点:float, 複素数:complex といった数値基本型やリスト:list, タプル:tuple, 辞書:dictionary と言った集まりを表現する基本型があります。

これらの基本型変数を処理する if then else 構文や loop 構文は C や Java に近いものです。

関数呼び出しについて C や Java と微妙に異なっています。

duck typing

Python の基本的な考え方として duck typing と言われる考え方があります。「があがあ」と鳴き「よちよち」と歩くのがアヒルならば、「があがあ」と鳴くメソッド「よちよち」と歩くメソッドを備えたものはアヒルと思ってしまえと言われます。C++/STL での考え方をさらに広げています。

solvePDE(..) 内部で使われる行列は、要素をインデックスでアクセスでき、それらの要素に係数を掛け、また加減算しているだけです。

Duck typing の考え方を solvePDE(..) 関数に適用しています。要素をインデックスでアクセスでき、それらの要素に係数を掛けることができ、また加減算できれさえすれば、それが何であっても同じ solvePDE(..) python プログラムを使えるように記述してあります。

Duck typing の考え方を適用することで、一つの solvePDE(..) 関数だけで大部分の偏微分方程式を扱えるように記述してあります。Python の凄さ、記述しやすさの一端を理解してもらえると思います。

とりあえず solvePDE(..) を使って偏微分方程式を解くだけならば duck typing の考え方まで理解する必要はありません。solvePDE(..) が任意の変数の偏微分方程式であっても、その値が任意次元のベクトル/テンソルであっても solvePDE(..) が使えることだけを覚えておいてください。

コメント

Perl などと同様に # 以降の文字列はコメントを意味します。

//@@
print "Hello World" # 定番のコードです
//@@@

bool int, float, complex, string 型

Python における bool, int, float, complex, string 型変数の注意すべき点を書いておきます。

bool 型

bool 変数は True または False のいずれかの値を取ります。Python では bool 値を返す演算子として == 演算子と is 演算子の二種類があります。

== 演算子で比較するとき、 0 や 0.0 と False は equal です。 1 や 1.0 と True は equal ですが 0.5 とは equal でありません。次のような具合です。

//@@
print 0 == False
print 0.0 == False
print 1   == True
print 1.0 == True
print 0.5 == True
//@@@
True
True
True
True
False

is 演算子は == 演算子より厳しく比較します。次のような具合です。

//@@
trueAt = True
falseAt = False
print 0 is False
print 0.0 is False
print 1   is True
print 1.0 is True
print 0.5 is True
print trueAt is True
print falseAt is False
//@@@
False
False
False
False
False
True
True

偏微分方程式を解く漸化式:fncRecurring(..) で指定されたインデックスが計算すべき領域であるか否かは、boundaryValue 引数値が True であることで判断します。このとき is 演算子を使います。== 演算子を使えません。boundaryValue が 1 の整数値または 1.0 の浮動小数点値かもしれないからです。

int 型

C でも同じですが int どうしの割り算は int に丸められます。// 演算子を使うことで丸め処理を明示できます。

//@@
print 5/4
print 5//4
print 5/4.0
print 5.0//4.0
//@@@
1
1
1.25
1.0

python 2.4 からは整数は任意サイズの big number も含むようになりました。

//@@
print 2**100
print 2**100/3.0
//@@@
1267650600228229401496703205376
4.22550200076e+029
float 型

double サイズの浮動小数点で扱われること以外に特に説明することはありません。

complex 型

デフォルトで浮動小数点の複素数:complex を扱えます。数値の直後に j を付けて複素数値を表現します。数値の付かない j だけでは複素数にできません。整数の複素数はありません。

//@@
print 2+3j
print 2+0j
#print 2+j  #NameError: name 'j' is not defined
print (2+3j)/2
//@@@
(2+3j)
(2+0j)
(1+1.5j)
string 型

「" .... "」 と " で囲まれた文字列が string 型です。「' .... '」と ' で囲まれた文字列も string 型です。文字列の途中で " が入るときなどでの記述を容易にするため二種類の文字列表記を可能にしています。「""" .... """」 と """ で囲まれた文字列も string 型です。このときは \n マークなしで改行を挿入できます。

//@@
print "abc"
print 'ab"ABC"c'
print """abc
    XYZ"""
//@@@
abc
ab"ABC"c
abc
    XYZ
漢字の扱い

Python プログラムの中で日本語文字列を使いたいとき、Shift JIS 日本語のコメント文字を使いたいときは、二行目までに「# coding=shift_jis"」の行を挿入します。

//@@
# 07.10.10
# coding=shift_jis
print "日本語を出力" # 日本語のコメント
//@@@

リスト:list, タプル:tuple, 辞書:dictionary

Python でのオブジェクトの集まりを表す「辞書:dictionary」「リスト:list」「タプル:tuple」といった基本型は、C や Java ではコンテナ・ライブラリとして実装されているものに似ています。基本型に組み込んでいるだけあって、これらは C/Java におけるコンテナ・ライブラリよりも使いやすく機能も豊富です。

リスト:list

Python におけるリストは C 言語における配列に近いものです。下のような C 言語におけるポインタで繋いだリストではありません。

struct StList{
    int m_element;
    struct StList* m_pStNext;
};

Python におけるリストは [ と ]で囲んであらわします。下のように使います。

//@@
# coding=shift_jis
print [1,2,3]
lstAt = [1,2,3]+[4,5,6]     # リスト同士を繋ぎます
print lstAt
lstAt.insert(3,100)         # インデックス 3 の位置に数値 100 の要素を挿入します。
print lstAt
print lstAt[3]  # index 3 の要素を取り出します
print lstAt[3:] # index 3 以降の全ての要素を取り出します
print lstAt[-1] # lstAt の最後の要素を取り出します
//@@@
[1, 2, 3]
[1, 2, 3, 4, 5, 6]
[1, 2, 3, 100, 4, 5, 6]
100
[100, 4, 5, 6]
6

そのほかにも python のリストで数多くの知っておくべきことがあるのですが、solvePDE(..) を使う分に必要なのは lstAt[一個の index] と要素をアクセスできることだけです。

タプル:tuple

Python におけるタプルは要素を変更できないこと以外は list に良く似ています。 ( と ) で囲んであらわします。

//@@
# coding=shift_jis
print (1,2,3)
tplAt = (1,2,3)+(4,5,6)
print tplAt
#tplAt.insert(3,100)
print tplAt
print tplAt[3]
print tplAt[3:] # index 3 以降の全て
print tplAt[-1] # lstAt の最後
print hash(tplAt)   # タプルの hash 値
//@@@
(1, 2, 3)
(1, 2, 3, 4, 5, 6)
(1, 2, 3, 4, 5, 6)
4
(4, 5, 6)
6
-319527650

python のタプルでも数多くの知っておくべきことがあるのですが、solvePDE(..) を使う分に必要なのは tplAt[一個の index] と要素をアクセスできることと、次に説明する辞書でタプルをキーに使えることだけです。リストは辞書のキーに使えません。

なお python におけるリストとタプルが よく似ていることも事実です。どちらで書いてもいいときが良くあります。このタプル/リストどちらでもよいときにはシーケンスという言葉で説明することがよくあります。このことも覚えておいてください。

辞書:dictionary

辞書型変数はキーと値の組にした要素の集合です。Perl での連想配列、C++/STL における map のようなものです。任意の python オブジェクトを値につかえます。

辞書は { と } で囲んだ間に key:value の組を「,」デリミタで連ねて表現します。要素へのアクセスは辞書変数ラベルの次に [key] を続けることで行います。下のように使います。

//@@
dctAt = { "test":100, (1,2):1000}
dctAt[1]="ABC"
dctAt['a']="abc"
dctAt[(1,2,3)] = 10**3
print dctAt['a']
print dctAt[(1,2,3)]
print dctAt
//@@@
abc
1000
{'test': 100, (1, 2): 1000, (1, 2, 3): 1000, 'a': 'abc', 1: 'ABC'}

辞書のキーにタプルを使うとき、dctAt[(1,2,3)] を dctAt[1,2,3] と記述することもできます。Python が旨く作られていると感じるところです。

この辞書のキーを行列やテンソルのインデックスのようにも記述できることを solvePDE(..) では有効に活用しています。辞書のキーを行列のインデックスのように記述できることは確りと覚えておいてください。下のような具合です。

//@@
dctAt = {}
dctAt[0,0]=1
dctAt[0,1]=2
dctAt[0,2]=3

dctAt[1,0]=4
dctAt[1,1]=5
dctAt[1,2]=6


dctAt[2,0]=7
dctAt[2,1]=8
dctAt[2,2]=9

dctAt[1,1]=100
print dctAt
//@@@
{(0, 1): 2, (1, 2): 6, (0, 0): 1, (2, 1): 8,
 (1, 1): 100, (2, 0): 7, (2, 2): 9, (1, 0): 4, (0, 2): 3}

モジュール

Python のモジュールは C 言語でのヘッダ・ファイルとライブラリの両方の性質を兼ね備えたものです。Python モジュールを import することで、そのモジュールに実装してある関数やクラスなどを利用できるようになります。solvePDE(..) 関数は sf.py モジュールに実装されています。

モジュールが大きくなって一つのファイルには収まりきらなくなり複数のファイルをディレクトリに収めるようになるとパッケージと言われるようになります。sf.py は solvePDE(..) 関数を含むモジュールですが、3D 表示機能をまとめた vpython はパッケージで実装されています。

モジュールやパッケージのインストール

Python でモジュールをインストールすると python があるディレクトリの下の Lib ディレクトリ以下にパッケージやモジュールが入れられます。通常はインストーラーが Lib 以下の配置場所を判断してインストールします。

solvePDE(..) 関数を含む sf.py モジュールはインストールするような汎用的なものではないので、偏微分方程式を解く作業をするカレント・ディレクトリに置いておきます。

モジュールの import

既存のモジュールやパッケージは import 文で指定することで他のモジュールから利用できるようになります。「from moduleName import *」とする方法と「import moduleName as aliasName」とする方法があります。名前空間を汚さないために「import moduleName as aliasName」と書くことを強く勧めます。下のような具合です。

//@@
# coding=shift_jis
import visual as vs
print vs.cross( [1,2,3], [4,5,6])   # ベクトルの外積を計算させます

import sf
#まだ dctBoundary などの引数にするオブジェクトがないのでコメント・アウトしておきます。
#sf.solvePDE(dctBoundary, fncRecurring, inIteration)
//@@@
<-3, 6, -3>

sf のように短いモジュール名のときは alias 名を設定せずに使います。タイプの煩わしさを厭わなければ as 以下の visual モジュールでも alias 名を使わないことも可能です。

モジュール実行

python のモジュールを実行するには、python プログラムが書きこまれた、拡張子が py のファイルを作り、それを python.exe の引数に指定します。次のような具合です。

//@@
print "Hello World"
//@@@
//copy \#####.### temp.py /y

「print "Hello Wold"」が書かれた temp.py ファイルをカレント・ディレクトリに用意できたら、それを下のように python.exe で実行します。

python.exe temp.py
Hello World

# .exe の拡張子名は省略することもできます。
python temp.py
Hello World

Python をインストールしたときインストーラーが python.exe のあるディレクトリに path を通しているので、python.exe ファイルのディレクトリ位置を指定する必要はありません。

バージョン違いなど複数の python をインストールしているときは full path 名で python.exe を指定します。

\lng\python25\python temp.py
Hello World

temp.py がカレント・ディレクトリ以外にあるときは temp.py もフル・パス名を使って表示せねばなりません。

python \my\vc7\mtCm\temp.py
Hello World

ただし python モジュールが python library 内にあるときは -m オプションを付けることで full path 名にしなくても、すむようになります。このときは .py 拡張子名も省略できます。

python -m pydoc bisect.bisect
Help on built-in function bisect in bisect:

上で実行した pydoc.py は python の標準ライブラリに入っているモジュールです。次に説明する introspection 機能を実現します。

pydoc.py モジュール, 内省:introspection

Python では 内省:introspection などと高級そうな言い方をされますが、ようはヘルプ文字列の表示です。Python 関数宣言の直後に書かれた """ ..... """ 文字列を表示する機能です。下のように testF() 関数にコメント文字が入っていたとき、これを pydoc.py モジュールで取り出せます。

//@@
def testF():
    """ test comment at testF() function
    """
    pass

def test2F():
    """ Now this is test2 comment at testF2()
    """
    pass
//@@@
//copy \#####.### temp.py /y

この temp モジュールの testF 関数のコメントが下のように見えます。

python -m pydoc temp.testF

Help on function testF in temp:

temp.testF = testF()
    test comment at testF() function

下のように temp モジュールだけを指定したときは、モジュール内全ての関数のコメントが下のように表示されます。

python -m pydoc temp
Help on module temp:

NAME
    temp

FILE
    c:\my\vc7\mtcm\temp.py

FUNCTIONS
    test2F()
        Now this is test2 comment at testF2()
    
    testF()
        test comment at testF() function
math モジュールの introspection

一般のモジュールは python ライブラリ Lib 以下のどこにあるか分りません。フルパスでのモジュール・ファイル名を pytdoc.py 引数に与えれば instrospect できます。でもカレント・ディレクトリまたは Lib 以下にあるモジュールについてはフルパスをかかなくてもファイル名だけで pydoc が探し出してくれます。偏微分方程式の数値解でよく使う math モジュールについて下のように introspect できます。

python -m pydoc math
Help on built-in module math:

NAME
    math

FILE
    (built-in)

DESCRIPTION
    This module is always available.  It provides access to the
    mathematical functions defined by the C standard.

FUNCTIONS
    acos(...)
        acos(x)
        
        Return the arc cosine (measured in radians) of x.
    
    asin(...)
        asin(x)
        
        Return the arc sine (measured in radians) of x.
        ・
        ・
    tan(...)
        tan(x)
        
        Return the tangent of x (measured in radians).
    
    tanh(...)
        tanh(x)
        
        Return the hyperbolic tangent of x.

DATA
    e = 2.7182818284590451
    pi = 3.1415926535897931

Introspection 機能はマニュアルを開かずにキーボードで数文字をたたくだけで help を表示してくれ便利です。キーボードからマウスに持ち帰る必要さえありません。ただ python -m pydoc を何度も打つのが面倒なので私は下のような bat ファイルを作って path の通ったディレクトリに置いています。

type \utl\pd.bat
python -m pydoc %1

下のような具合に使っています。

pd math.log

Help on built-in function log in math:

math.log = log(...)
    log(x[, base]) -> the logarithm of x to the given base.
    If the base not specified, returns the natural logarithm (base e) of x.
pdb.py モジュール, デバッガ

後回し

if 構文

if else 構文は C 言語と殆ど同じです。Curry Brace 「{ と }」が使われないことと 「else if」の代わりに elif が使われることが C 言語と違います。

Curry Brace 「{ と }」が使われない代償として if 条件文の最後に ":" が必要なことと、インデントが強制される点が C 言語と大きく違う点でしょう。次のような具合に記述します。

//@@
inAt = 3

if inAt == 0:
    print "0"
elif inAt == 1:
    print "1"
else:
    print "others"
//@@@
others

インデントの強制は Python の特徴の一つです。C 言語のように複数のインデント・スタイルが共存していないため、他人のソースを読むとき楽をできます。

なお、python では C 言語での switch 文が存在しません。上のように if ... elif ... elif ... else ... 文で記述します。

loop 構文

C 言語での下の for loop 構文を考えて見ましょう。

for (int i=0; i<10; ++i){
    cout << i << endl
}

上の C コードと同じ内容を python では次のように書きます。

//@@
for i in range(0,10,1):
    print i
//@@@

python temp.py
0
1
2
3
4
5
6
7
8
9
range(..) 関数は introspection により次のようなものであると分ります。
pd range

C:\my\vc7\mtCm>python -m pydoc range 
Help on built-in function range in module __builtin__:

range(...)
    range([start,] stop[, step]) -> list of integers
    
    Return a list containing an arithmetic progression of integers.
    range(i, j) returns [i, i+1, i+2, ..., j-1]; start (!) defaults to 0.
    When step is given, it specifies the increment (or decrement).
    For example, range(4) returns [0, 1, 2, 3].  The end point is omitted!
    These are exactly the valid indices for a list of 4 elements.

代入文

C 言語とは異なり代入文は値を持ちません。if 文の条件式に代入文を書いても、下のようにエラーになるだけです。これにより C 言語でよくやる間違いを犯さずにすみます。

//@@
inAt = 3
if inAt = 5:
    print "5"
//@@@

python temp.py
  File "temp.py", line 2
    if inAt = 5:
            ^
SyntaxError: invalid syntax

C 言語とは異なり代入文の左値に複数の変数を記述できます。これはプログラムの記述を簡便にしてくれます。下のような具合です。

//@@
x=3
y=5
x, y = y, x # exchange
print (x,y)

i, j = (10,11) # multiple assignment
print (i,j)
//@@@

python temp.py
(5, 3)
(10, 11)

sf.solvePDE(..) が呼びだす fncRecurring は下のような偏微分方程式の漸化式を処理する関数です。ここで index は tuple で渡されます。偏微分方程式が三変数のときは index の tuple の長さは 3 になります。二変数のときは index の tuple の長さは 2 になります。実際の対象に合わせて、下のように index を個別変数 j,k に割り振りなおします。

def fncPDE(dctRslt, index, boundaryValue):
    j,k=index
    dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                          + dctRslt[j+1,k]
                          + dctRslt[j,k-1]
                          + dctRslt[j,k+1])

上のように代入文の左値に複数の変数を置けることを利用して一つの sf.solvePDE(..) 関数だけで、任意個数の変数の偏微分方程式を扱えるようにしています。

デフォルト引数とキーワード引数

C++ 関数でのデフォルト引数値は右側の連続する変数値をまとめてデフォルト値にするものです。引数の順序を変えられないので、デフォルト値を選択できません。

int testF(inLeftAg=1, inRightAg=4)
{
    return inLeftAg+inRightAg
}

cout << testF(1) << endl;   // inRightAg==4 であり 5 をプリントする
cout << testF(1, 7) << endl;   // inLeftAg のデフォルト値を使えない。

でも python では引数変数名を明示的に記述することで、残りのデフォルト値を有効にできます。キーワード引数名を指定することで、デフォルト以外をとる引数を指定できます。下のような具合です。

//@@
def testF(inLeftAg=1, inRightAg=4):
    return inLeftAg+inRightAg

print testF(3)  # same as testF(3,4)
#print testF(?, 5)  # can't call at C
print testF(inRightAg=5)    # same as testF(1, 5)
//@@@
7
6

solvePDE(..) 関数は下のデフォルト引数値を設定しています。

def solvePDE(dctBoundary, fncRecurring, inIteration=64, dfltValue = 0, arRslt=None):

dfltValue 引数名を指定することで、他のデフォルト値を有効にしたまま、dfltValue 引数値をベクトル配列 sf.Sc.array([0,0,0]) に変更するようなことが可能です。これは div E = 0 のようなベクトル値を持つ偏微分方程式で必要となります。

        ・
solvePDE(dctBoundary, fncRecurring, dfltValue = sf.Sc.array([0,0,0]))
        ・

変数の生成消滅の管理はレファランス・カウント

Python の全てのオブジェクトの生成消滅はレファランス・カウント方式で行われます。C プログラマーには違和感があるでしょうが、下のコード例のように関数の内部でオート変数として作られた辞書変数:dctRslt を戻り値にできます。dctRslt オブジェクトが delete されるのは dctRslt が他から一切参照されなくなったときです。

def solvePDE(dctBoundary, fncRecurring, inIteration=64, dfltValue = 0, arRslt=None):
    ・
    ・
    dctRslt={}
    ・
    ・
    return dctRslt

さらに詳しく python を勉強してみる

後回し

vpython 文法:solvePED(..)を扱うために必要な最小限の文法解説

とりあえずは、vpython について何も知らなくても sf.make3dRowCol(dctRslt) を呼び出すことさえ知っていれば計算結果を表示できます。ここでの vpython に関しての必要な知識は 3D 表示された図形の操作方法とインストールについてだけです。

vpython 3D 図形の回転と拡大

vpython が表示する 3D 図形の左右回転は、マウスの右ボタンを押したまま左右にドラッグすることで行います。 3D 図形の上下回転は、マウスの右ボタンを押したまま上下にドラッグすることで行います。3D 図形の拡大・縮小はマウスのセンター釦を押したま上下にドラッグすることで行います。

vpython のインストール

Enthought パッケージをインストールして python2.4 の方は ここの VPython-Win-Py2.4-3.2.9.exe をダウンロードして実行してください。

Python2.5 向けの vpython もありますが機能に大差ありません。マニュアルは 2.4 向けと同じです。

sf モジュールを動かす分には python2.5 向けの vpython でも動きます。

sf 文法:ここでの偏微分方程式を扱うために必要な最小限の解説

sf は私の開発した行列・ベクトル計算を行わせるソフトです。エディタで書いている数式のままに行列ベクトル計算させます。コンソールから CUI で実行させます。計算結果は拡張子が val の sf 変数ファイルに残ります。この sf 変数ファイルは、sf 式の中の変数名として記述できます。

sf には slice 操作があり、行列テンソルの部分行列/ベクトルを抜き出せます。それらをグラフ表示、3D 表示させられます。膨大な数値の塊である偏微分方程式の数値解を具体的な形として人間がイメージするとき強力なツールとなります。特に 汎用的な可視化ツールが設けられない三変数以上の偏微分方程式の数値解のイメージをユーザーが作り上げるに特に便利なツールです。

sf.py には python の行列/テンソルを sf で計算扱えるファイル変数にする putSf(dctRslt), getSf("variableName") 関数を用意しています。これをコンピュータ用語を使って言い換えてみます、シリアライズ関数 putSf(.) getSf(.) を用意しています。これはファイルとの間の出し入れだけでなく、行列処理・視覚化も可能にするシリアライズです。

コンソールから下のように実行させます。

sf.exe "3+4"
< 7 >
sf "mTemp =`σx^3+ 5`σx^2 +`σx"
< 5, 2 >
< 2, 5 >
sf `σx
< 0, 1 >
< 1, 0 >

`σx は sf20.ini というテキスト・ファイルの中で宣言してあるパウリ行列です。

mTemp = ... と sf 式の左側の変数名については、カレント・ディレクトリに、拡張子が .val の同じ名前のファイルを作ります。(複数の sf 式があるときは最後の式についてだけ sf 変数ファイルを作るように、sf2.0 から改良しました。)。この変数ファイルは下のように、再度 sf 式の中で参照できます。カレント・ディレクトリに左値が残るので、試行錯誤しながら、いろんな試し計算をするのに便利です。

type mTemp.val
% BaseDoubleMatrix 2, 2
< 5, 2 >
< 2, 5 >

sf "vTemp =<1,3+4i>, mTemp vTemp"
<  11 +8i, 17 +20i >

左値がない sf 式では _dt.val の sf ファイル変数がカレント・ディレクトリに作られます。

type _dt.val
% BaseComplexDouble 2
< 11 +8i, 17 +20i >

この _dt.sf 変数ファイルも次の sf 式の対象にできます。下のような具合です。

sf "_dt*2"
< 22 +16i, 34 +40i >

sf のインストールについて

sf 自体は有償のソフトですが、ここから評価版をダウンロードできます。評価版といっても機能制限はありません。少し遅くなるのと、権利表示が出てくるのがうざいのと、Vc7VrfyMDRt110D.zip ファイルがキー・ファイルとして必要な点で異なっているだけです。

なお sf20beta.zip の側をダウンロードしてください。少し問題が残っていますが、ここでの偏微分方程式を扱う程度の使い方では問題ありません。

sf の help

/? を引数に sf を実行させることで sf の使い方を表示させられます。下のような具合です。

sf /?
Your machine code is : 4b012029060cb1ed
sf version 2.0β.  Build::Vc7:Vc1300.9466 vrfyPrML10C_298947b17a.lib 2005.05.10

sf.exe::06.05.15 build
::four operations: + - * / ^ >sf "3+4+5, 3-4, 3*(4+5), (3+4)/5,3^0.4"
    default calculation result is in _dt.val. Use the last result. >sf _dt+1+2
::complex number operations: >sf "3+4, 3i-4,3*(1+4i),(4+i)/3,i^0.5+i"
::transform to hex integer   >sf 0x=1024+0x10+0b1111_0101
::vector/matrix expression   >sf "vector=<0,1,2,3,4>,vector[0]=1000"
  generate size 10, 0 element vector  >sf "vector=<<10>>"
  generate Matrix of column/row e 3/4 >sf "matrix=[[3,4]], matrix[0,1]=1000"
  generate 3x3 square matrix & set row >sf "matrix=[[3]], matrix[0,*]=<1,2,3>"
  matrix row and diagonal access      >sf "matrix[*,2]=<4,5,6>,matrix[*,*]=<7,8,9>"
:: inner product             > sf " >, "

::using slice array -- modification by 3 elements vector:
    vector 0,2,4 element  >sf "vector=<<9>>,vector[<0,3,2>] = <1,2,3>"
                  result  < 1, 0, 2, 0, 3, 0, 0, 0, 0 >
    matrix diagonal element by slice >sf "matrix=[[3]], matrix[<0,3,4>]=<1,2,3>"
        result  < 1, 0, 0 >
                < 0, 2, 0 >
                < 0, 0, 3 >
::repeat with slice parameter and generate vector
    generate vector by temporary parameter  >sf "<<0,4,1 @t | t^2>>"
         result < 0, 1, 4, 9 >
::options         
    /s  : silent    /ns: not silent at console
    /p N: print digit number
    /zs float int:  zero supress, vector size
    /pvt float   :  tiny pivot value to invert matrix

::embedded functions         
        example     >sf !log(1+i)
        result      < 0.346574+0.785398i >
  basic functions
    !sin  !cos  !tan  !sinh !cosh !tanh !exp !log !asin !acos !atan !log10
  other functions
    !abs:   absolute for each element
    !norm:  square root norm
    !fft:   Fast Fourier Transformation
    !rft:   Reverse Fast Fourier Transformation
    !dggr:  dagger for matrix, conjugate for vector and scalar
    !image, !real:  get image/real from complex number
    !lt:    less than: bool lt(left,right)
    !shift: shift vector/matrix data:shift(data, inSift), inSift>0 ==> right sift
    !size:  get vector size
    !sqrt : square root
    !sum  : sum up all vector element
    !prdct: product of all vector element

sf のループ構文

行列・ベクトル計算を行わせる sf.exe には if 分がありません。行列・ベクトル計算では if 文を必要とする複雑さにならないことが大部分だからです。どうしても if 文が必要なときは python script を呼び出すことでできるので、そちら側で if 文を含む複雑な処理を行います。

でも繰り返しは頻繁に必要になります。そのため下のような形式での loop 構文を使えるようになります。

<>

k ≡ start + n * stride  但し n ∈[0,size)

上で k には start + n * stride で決まる長さが size の等差数列の数値が入ります。size は整数値でなければなりませんか、start, stride は同じサイズのベクトルや行列も使えます。<<... @k| sf 式 >> の結果は sf 式の値が行列のときはベクトルになります。sf 式の値がベクトルのときは行列になります。下のような具合に sin 関数のベクトルを作れます。

sf "N=64, <<0,N,1@k|!sin(2`π k/N)>>"
<            0,    0.0980171,      0.19509,     0.290285,     0.382683,     0.471397,
       0.55557,     0.634393,     0.707107,      0.77301,      0.83147,     0.881921,
       0.92388,      0.95694,     0.980785,     0.995185,            1,     0.995185,
      0.980785,      0.95694,      0.92388,     0.881921,      0.83147,      0.77301,
      0.707107,     0.634393,      0.55557,     0.471397,     0.382683,     0.290285,
       0.19509,    0.0980171, 1.22461e-016,   -0.0980171,     -0.19509,    -0.290285,
     -0.382683,    -0.471397,     -0.55557,    -0.634393,    -0.707107,     -0.77301,
      -0.83147,    -0.881921,     -0.92388,     -0.95694,    -0.980785,    -0.995185,
            -1,    -0.995185,    -0.980785,     -0.95694,     -0.92388,    -0.881921,
      -0.83147,     -0.77301,    -0.707107,    -0.634393,     -0.55557,    -0.471397,
     -0.382683,    -0.290285,     -0.19509,   -0.0980171 > 

計算結果のグラフ表示

計算結果は sfGraph.py に渡すことでグラフにできます。sf 変数ファイルが _dt.val のときは、その変数名を省略できます。なお Windows では Enthought python をインストールすると拡張子 py が python.exe と関連付けられます。sgGraph.py をコンソールから打ち込むと python.exe sgGraph.py が実行されます。

sfGraph.py _dt
または
sfGraph.py

sfGraph.py を一々打ち込むのが面倒なので、私は下の pg.bat を使っています。

type pg.bat
sfGraph.py %1

sf は私自身が欲しくて作ったソフトです。自分では一番長い時間使っているソフトです。計算結果が HDD に残るので、プログラムを組み必要がありません。計算結果に応じて別の計算をさせる試行錯誤の連続の過程で使っています。皆様にもお役にたてると思います。登録料 5000 円の有償ですが、電卓代金だと思ってください。

sf.exe のマニュアル

sf のより詳細なマニュアルがここにあります。sf ver1.2 向けなのですが、sf ver2.0 は殆ど upper compatible にしてあるので入門的な使い方をするぶんは十分です。


sf.py モジュールについて

sf.py モジュールには sovePDE(..) 関数も含め、数値計算で頻繁に使う機能をまとめてあります。そのうち偏微分方程式の数値解を求めるときに必要になる関数について説明します。

sf.py のインストール

sf.py は特に偏微分方程式専用のモジュールであり python 自体にはインストールしません。この sf.py をカレント・ディレクトリに置いてください。

sf.mrng

Python の range 組み込み関数を使って多重ループを何度も記述するのが面倒なので mrng 関数を作りました。

//@@
dctAt = {}
for i in range(3):
    for j in range(3):
        dctAt[i,j] = i^2 + j^2
print dctAt
//@@@
{(0, 1): 1, (1, 2): 7, (0, 0): 0, (2, 1): 3, (1, 1): 0,
 (2, 0): 2, (2, 2): 4, (1, 0): 1, (0, 2): 6}

上の二重ループを python code を下のように一回のループだけで書けます。sf.mrng はリカーシブに記述してあり、何重ループでも同様に記述可能です。何度もループ構文を扱う偏微分方程式の数値解プログラムで mrng を導入するのも理解してもらえると思います。

//@@
import sf
dctAt = {}
for i,j in sf.mrng(3,3):
        dctAt[i,j] = i^2 + j^2
print dctAt
//@@@
{(0, 1): 1, (1, 2): 7, (0, 0): 0, (2, 1): 3, (1, 1): 0,
 (2, 0): 2, (2, 2): 4, (1, 0): 1, (0, 2): 6}

ループ構文のような基本的な記述さえ、ユーザーの望むように改修できるのが python のすばらしさです。mrng のソースも読んでもらうことを希望します。Python の素晴らしさを味わってもらえるはずです。

sf.arSqnc

浮動小数点、複素数、ベクトル/行列の等差数列を表すのに使います。Python の range(.) 関数では整数の等差数列しか表現できません。range(start=0, stop, stride=1) の引数を渡すとき、大小関係を比較する stop 引数を渡しているので複素数の等差数列さえ表現できません。

次のように使います

//@@
import sf
lstAt=[]
for val in sf.arSqnc(1+2j, 5, 1j):
    lstAt.append(val)
print lstAt
//@@@
[(1+2j), (1+3j), (1+4j), (1+5j), (1+6j)]

下のような説明をソースにつけてあります。

pd sf.arSqnc

sf.arSqnc = arSqnc(startOrSizeAg, sizeAg=None, strideAg=1)
    等差数列 を start, size, stride 引数で生成する
    arithmetic sequence is generated with argment start, size, stride

sf.masq

arSqnc の多次元版です。

sf.make3dRowCol

行列データより vpython を使って 3D 図形を表示します。ただし行列データを左側に 90 度回転させた状態で表示増す。

行列の column, row インデックスと座標位置ベクトル (x,y) とは配置順序が逆です。旨く対応しません。行列の座標データを紙の上に一覧表で書くとき、[0,0] インデックスは左上に書かれます。でも (0,0) 座標は左下にかかれます。

     y
     │Matrix[0,0]
     │  ----- row -------- 
     │ |                   
     │ |Matrix(column,row) 
     │column               
     │ |                   
     │ |  posVector = (x,y)
   ─┼───────── x 
     │pos(0,0)

この分りにくさの対策として、make3dRowCol(.) に渡された行列データは、下のように 90 度左側に回転させて表示します。

     y
     │
     │ |      x
     │ |      i
     │row     r
     │ |      t
   Matrix[0,0] M
     │  --- column ----
   ─┼───────── x
     │

その他の機能として、下のようなオプション表示設定がありますが、とりあえずは これらを弄らなくても偏微分方程式の行列データを表示できます。

make3dRowCol(mtrxAg, blMesh = True, blAxis = True, colourRgb = vs.color.blue
        , blDisplayMinMaxLabel = False, xyRatio = 1.0):

あと make3dRowCol(..) を走らせると mtrxAg 引数の最大値、最小値と対応する column/row index をコンソールに下のように表示します。上の blDisplayMinMaxLabel キーワード引数を True にしてやれば 3D vpython 図形の中にも min/max 値を表示します。でも、デフォルトで表示するコンソール値の方が便利でしょう。

index:(19, 14)  max:     1
index:(29, 15)  min:-0.0532405

sf.putSf

偏微分方程式の計算データ行列/テンソルを sf ファイル変数としてカレントディレクトリに書き出す関数です。デバッグ途中に便利に使えます。デフォルト sf ファイル変数名は _rs です。

def putSf(valAg, fileNameAg="_rs"):    

sf.getSf

後回し

偏微分方程式の数値解法

差分法アルゴリズム

ここでの偏微分方程式で使うアルゴリズムは、行列の Gauss Seidel 法のように漸化式を何度も適用して偏微分方程式を満たすように数値を修正していくものです。偏微分方程式に対応する漸化式の作り方が肝となります。

どのように偏微分方程式に対応する漸化式を作るかは、私の下手な説明より、ここにあるラブラス方程式波動方程式の説明を呼んでもらったほうが分ると思います。

ここの説明を読んでもらえば、ラブラス方程式、波動方程式以外の偏微分方程式に対応する漸化式も自分で作れると思います。 sf.solvePDE(.) python 関数

solveLaplacian(..), solvePDE(..), itrSlvPDE(..) python 関数

今までの知識を前提に、やっと solveLaplacian(..), solvePDE(..), itrSlvPDE(..) 関数の使い方の厳密・詳細な説明に入れます。

境界条件行列

sf.py モジュールを使って偏微分方程式を解くためには境界条件を表す辞書行列(以下では境界条件辞書行列と呼びます)を作らねばなりません。それを solveLaplacian(..), solvePDE(..), itrSlvPDE(..) python 関数 の dctBoundary 引数に与えることで漸化式の繰り返し処理が行われます。

境界条件辞書行列は sf.py モジュールが扱える偏微分方程式の定義域も定めます。その定義域は n 次元直方体領域に限られます。n 次元直方体領域を同じサイズの小さな n 次元直方体で分割し、そのメッシュ点上の値を数値解として計算します。円形領域は扱えません。二変数のときは長方形の領域のメッシュ上の数値解を求めることになります。

境界条件辞書行列:ここでは より具体的に solvLapalacian(..), solvePDE(..), itrSlvPDE(..) の dctBoundary 行列引数の行列サイズは、メッシュ点上の数値解として求めようとしている行列と同じサイズの辞書行列となります。(逆に言えば dctBoundary の行列 column/row サイズより、偏微分方程式の数値解を保存する行列 column/row などのサイズを sf.py 側で定めます。) dctBoundary 行列は各メッシュ点の性質を表現する辞書行列だともいえます。dctBoundary 辞書行列の各要素には True, False, 文字または境界値のいずれかの値を設定します。

最初の「 solveLaplacian(..) の使い方の概要」で使った境界条件行列に下のように False を追加してみましょう。

True True True True True True True True True True  
                                                   
True True True True True True True True True True  
                                            
True True True True True True True True True True  
                                                   
True True True 2    2    2     2   True True True  
                                                   
True True True 2   False False 2   True True True  
                                                   
True True True 2   False False 2   True True True  
                                                   
True True True 2    2    2     2   True True True  
                                                   
True True True True True True True True True True  
                                                   
True True True True True True True True True True  
                                                   
True True True True True True True True True True  

False が設定されたメッシュ点はデフォルト値 0 に設定されます。プログラムと結果は下のようになります。

//@@
import sf
N = 30
N =  10
#N =  8
dctBoundary={}
for index in sf.mrng(N,N):
    dctBoundary[index]= True

for i in range(N/2-2, N/2+2):
    dctBoundary[N/2-2,i]=2;
    dctBoundary[N/2+1,i]=2;
for i in range(N/2-1, N/2+1+1):
    dctBoundary[i,N/2-2]=2;
    dctBoundary[i,N/2+1]=2;

for index in sf.masq([N/2-1,2], [N/2-1,2]):
    dctBoundary[index]=False;


# solve Partial Differential Eqation
dctRslt = sf.solveLaplacian(dctBoundary)
sf.make3dRowCol(dctRslt)                                # visualize the result
sf.putSf(dctBoundary)
sf.putSf(dctRslt,"rslt")
//@@@

中心の False の領域がデフォルト値 0 に張り付いています。

なお、上のコードでは sf.mrng(..) や sf.masq(..) イタレータを使っています。これらを使って、二重ループの記述を一重の記述に簡略化させています。このイタレータの便利さを分ってもらえるでしょうか。

solveLaplacian(..) 関数

solveLaplacian(..) 関数について introspection を行います。

pd sf.solveLaplacian

python -m pydoc sf.solveLaplacian 
Help on function solveLaplacian in sf:

sf.solveLaplacian = solveLaplacian(dctBoundary, inIteration=64)
    'solve partial differential equation: solvePDE(dctBoundary, inIteration=64)
        
        dctBoudary argment: dictionary
            dctBoudary[i,j, .. , k] = True or False or Character
                True: calculate recurrence formula for Laplacian with dctRslt
                False:dctRslt[i,j, .. ,k] == 0, not calculating area
                Character: user extension <== don't use fo solveLaplacian
        inIteration argemnt: number of iteration times
    
        return: tensor dictionary with float value
    '
solveLaplacian(..) 関数は dctBoundary 引数以外に、漸化式をメッシュ点全体に計算させる繰り返し回数を設定できることが分ります。

sf.solvePDE(..) 関数

solvePDE について introspection を行います。

pd sf.solvePDE

python -m pydoc sf.solvePDE 
Help on function solvePDE in sf:

sf.solvePDE = solvePDE(dctBoundary, fncRecurring, inIteration=64, dfltValue=0, blInterpolation=True, arRslt=None)
    'solve partial differential equation by user defined recurring function.
        
        dctBoudary argment: dictionary
            dctBoudary[i,j, .. , k] = True or False or Initial value
                True: call fncRecurring function which will be called for user
                      with dctRslt and index argment
                False:dctRslt[i,j, .. ,k] == dfltValue, not calculating area
        
        fncRecurring: function defined by user
            argment fncRecurring(dctRslt, index)
    
        inIteration: number of iteration times
        dfltValue: defualt value at dctBoundary[index] == True or False 
        blInterpolation: True: enable linear interpolation at mesh edge
        arRslt: Use user defined mesh matrix/tensor
        return: matrix or tensor dictionary with float value
    '

上を見て分るように solvePDE(..) では境界条件辞書行列:dctBoundary の他に fncRecurring が必須の引数になっています。fncRecurrring 引数には偏微分方程式に対応する漸化式関数を指定します。この関数が solvePDE(..) から call back されます。fncRecurring(.) の詳細については少し複雑なので節を変えて再度説明します。

inIteration 引数の意味は繰り返し回数であり、solveLaplacian(..) のときと同じ意味です。

dfltValue 引数は、境界条件辞書が True/False になっている箇所のメッシュ点の初期値を定めます。通常はデフォルトの 0 で済むと思います。ベクトル値は偏微分方程式を扱うときは、dfltValue 引数に明示的に 0 ベクトルなりのベクトル値を設定せねばなりません。そうしないと漸化式がベクトル値と 0 scalar 値との計算をすることになるからです。

arRslt は計算結果を保存する行列引数をユーザー側で用意したいときのために設けられている引数です。通常は使いません。

blInterpolation 引数

blInterpolation 引数も新たに加わっています。これは bool 値であり、境界条件辞書の端面での線形外挿補間を solvePDE(..) 側で行うことを意味します。blInterpolation 引数値を明示的に False にすると、端面での線形外挿補間を行わなくなります。ユーザー側:fncRecurring 側で端面を判定し、端面処理を記述せねばなりません。FDTD への適用例で、これを False にする必要があります

なお blInterpolation 引数が True のとき dctBoundary 参照引数の端面に対応する要素は True から None に変更されます。

端面での線形外挿補間をデフォルトで自動的に行わせていることが sf.py モジュールでの偏微分方程式の自慢したい点です。これにより偏微分方程式が楽に解けるようになります。偏微分方程式の端面での処理を考えなくても済むからです。blInterpolation 引数を False にしたときは、ユーザーが自分で端面の処理を記述せねばなりません。多くのばあい特別な漸化式処理が必要となります。f(j-1,k-1) などのインデックスの着いた漸化式を使って偏微分方程式の数値解を求めるのですが、特別な処理をしてやらないと j などのインデックスが端面で 0 の値となったとき、j-1 が定義域を超えてしまうからです。反対側の端面では j+1 のインデックスが定義域を超えます。これらの処理を毎回考えねばならないのは面倒です。たいていの場合は線形外挿補間で済みます。この自慢の点を理解してもらえるでしょうか。

fncRecurring 引数

solvePDE(..) 関数から call back される漸化式関数 fncRecurring(..) は下の二つの引数を持たねばなりません。
def fncPDE(dctRslt, index):
    j,k = index
    dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                          + dctRslt[j+1,k]
                          + dctRslt[j,k-1]
                          + dctRslt[j,k+1])

上で fcnPDE は solvePDE(..) から三つの引数を設定されて呼び出されます。

dctRslt 引数は solvePDE(..) 内部で使っている偏微分方程式の数値解を保持する辞書行列変数です。 dctRslt 引数は辞書を表す変数であり参照渡しです。値渡しではありません。fncPDE(..) の処理中に修正する dctRslt 要素は sovePDE(.) 関数の中で保持している偏微分方程式の数値解を保持する辞書行列と同じオブジェクト・インスタンスです。

index 引数には漸化式 fncPDE(..) 関数に処理してもらいたい辞書行列:dctRlst のインデックス・タプルが設定されています。通常は辞書の row,col 値 j, k に分けなおします。

ここで注目して欲しいのは solvePDE(..) 関数が扱うサイズの自由性です。sovePDE(..) 関数は変数が三階以上の偏微分方程式も扱えます。div E == 0 のようなベクトルの偏微分方程式でさえ扱えます。C や Java といった言語では、こんな書き方はできません。Duck typing の python だからできる芸当です。

sf.itrSlvPDE(.) python generator

itrSlvPDE について introspection を行います。

pd sf.itrSlvPDE
python -m pydoc sf.itrSlvPDE 
Help on function itrSlvPDE in sf:

sf.itrSlvPDE = itrSlvPDE(dctBoundary, inIteration=64, dfltValue=0, blInterpolation=True, arRslt=None)
    'generatro function iterating for solving partial differential equation: 
        dctBoudary argment: dictionary
            dctBoudary[i,j, .. , k] = True or False or Initial value
                True: yield with dctRslt and index argment
                False:dctRslt[i,j, .. ,k] == dfltValue, not calculating area
        inIteration: number of iteration times
        dfltValue: defualt value at dctBoundary[index] == True or False 
        blInterpolation: True: enable linear interpolation at mesh edge
        arRslt: Use user defined mesh matrix/tensor

        yield dctRslt, index, boundaryValue
    '

solvePDE(..) と少し違うだけです。solvePDE(..) から fncRecurring 引数をとり、関数呼び出しをするかわりにイタレーターで dctRslt と index を yield しているだけです。イタレータ・ループに中に漸化式を記述することになります。Python に慣れている方は、solvePDE(..) よりも itrSlvPDE(..) の方が記述しやすいと思います。

個別の偏微分方程式への適用例

△φ = 0 :電位分布

最初は二次元での電位分布です。solveLaplacian(..) で解けるのですが、あえて solvePDE(..) を使います。自前で導いた漸化式を使う例として分りやすいからです。

△φ = 0 の漸化式

∂^2ψ/∂x^2, ∂^2ψ/∂t の二階微分は次のように離散化できます
∂^2φ/∂x^2[j,k] ≒ (φ[j+1,k] + φ[j-1,k] - 2*φ[j,k])/Δx^2 ------------ (1)
∂^2φ/∂t^2[j,k] ≒ (φ[j,k+1] + φ[j,k-1] - 2*φ[j,k])/Δt^2 ------------ (2)
0 == ∂^2φ/∂x^2 + ∂^2φ/∂t == (1)式 + (2)式 の関係は、Δx = Δy = h として単純化し、次のように変形できます。
0 == (1)式 - (2)式 
  == (φ[j+1,k] + φ[j-1,k] - 2*φ[j,k])/h^2 - (φ[j,k+1] + φ[j,k-1] - 2*φ[j,k])/h^2
∴
φ[j,k] = 1/4(φ[j+1,k] + φ[j-1,k] + φ[j,k+1]+φ[j,k-1])

ようは φ[j,k] の周囲四つの平均をとるということです。なだらかな表面を作っていくということです。

△φ = 0 数値解

漸化式が求まったら、境界条件を与えて solvePDE(.) を適用するだけです。二箇所の閉じた境界での電位を与えることにしましょう。丸い境界条件は指定が面倒なので、四角い境界条件にします。

//@@
# coding=shift_jis
# 07.10.28 偏微分方程式 ラプラスの式 ≡ ∂^2φ/∂y^2 + ∂^2φ/∂x^2 == 0 電場
# 30 x 30  matrix position [15,10] と [15,20] の位置に 3 position 正方形を置き
# 0`V と 1`V の電位とする
# 外周 [:,0], [0,:], [:,N-1], [N-1,:] の位置は free:None とする
# <== free の位置の値は gradient より interpolate する
import sf
N = 30
dctBoundary = {}    # 境界条件を設定する辞書

# とりあえず 30x30 行列の全部を 計算可能:True にする
for index in sf.mrng(N,N):
    dctBoundary[index] = True

# x,y position [10, 15] にある 0`V 電位の設定
for index in sf.masq( (9,3), (14,3) ):
    dctBoundary[index] = -1

# x,y position [20,15] にある 1`V 電位の設定
for index in sf.masq( (19,3), (14,3) ):
    dctBoundary[index] = 1

iterationAt = 40
def fncPDE(dctRslt, index):
    j,k=index
    dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                          + dctRslt[j+1,k]
                          + dctRslt[j,k-1]
                          + dctRslt[j,k+1])

dctRslt = sf.solvePDE(dctBoundary, fncPDE, iterationAt)

sf.make3dRowCol(dctRslt)
sf.putSf(dctRslt, "rslt")
//@@@

偏微分方程式を解くプログラムといっても大部分は境界値の設定です。慣れるまでは初期条件辞書行列 dctBoundary インデックス指定に戸惑うと思います。辞書行列が 90 度左回転させた状態になるからです。(x,y) 位置指定を中心にイメージしてもらうと早く慣れると思います。

偏微分方程式を実際に解いているのは 7 行の関数: fncPDE(..) です。漸化式を計算しているだけです。

表示処理は sf.make3dRowCol(dctRslt) の一行だけです。この関数呼び出しだけで下のような 3D 表示を行わせられます。この 3D 図形は、ぜひとも上の python コードを動かして実際に vpython で表示させた図形を自分の手で動かし、様々の角度から御自分の目で眺めてください。



sf.make3dRowCol(dctRslt) を呼び出すことで、コンソールには下のような偏微分方程式の解の最大値、最小値が表示されています。この最大値の位置に 3D 図形では小さな赤い球を置いてあります。最小値の位置には小さな緑の球を置いてあります。これらによって偏微分方程式の解の形が明確にイメージできるはずです。

index:(19, 14)  max:     1
index:(9, 14)  min:    -1

この偏微分方程式の数値解は、二本の割り箸を束にしてゴム膜に突き刺したような感じです。調和関数の初等的なことは知っているつもりでしたが、こんな様子まではイメージできていませんでした。

この計算結果の 3D 図形を見ていると Lagrange や Euler は、これらの解の様子を頭の中でイメージできていたのだろうと思えるのですが、ここらに詳しい方、如何でしょうか?

△φ = 0 計算結果の誤差の検討

上の python コードの最後に sf.putSf(dctRslt, "rslt") の行があるので、偏微分方程式の計算結果が入っている行列 dctRslt を rslt の名前の sf 変数ファイルとしてカレント・ディレクトリに出力されています。これを使って偏微分方程式数値解の精度について少し考えて見ましょう。

y = 14 の横軸の様子を下のように取り出せます。

rslt[*,14]
<  0.0584448,          0, -0.0592695,  -0.123497,  -0.196908,  -0.284261,  -0.391707,
   -0.528811,  -0.714532,         -1,         -1,         -1,  -0.651732,  -0.405227,
   -0.200059, -0.0105229,   0.180318,   0.389354,   0.642262,          1,          1,
           1,   0.711852,   0.525385,    0.38744,   0.277658,   0.185333,   0.103376,
  0.0262763, -0.0508239 >

この数値の塊のベクトルを下のようにグラフ表示させられます。上の sf 式には左値がないので _dt.val に上のベクトル値が記録されています。pg.bat を実行させれば、コマンド引数がないので _dt.val を表示します

pg

下のようなグラフになります


左右の両端で +0.05,-0.05 の値になっています。でも電磁気学からの直感では左側はマイナス値であり、右側は + 値のはずです。このことから、この偏微分方程式の計算精度は 5% を切っています。目の子で 1 割り程度の精度だと推測されます。

この偏微分方程式は上下に対象になります。rslt[*,14] の反対側も見てみましょう。

rslt[*,16]
<    0.05912,          0, -0.0599076,  -0.124662,  -0.198429,  -0.285925,  -0.393278,
   -0.530055,  -0.715237,         -1,         -1,         -1,  -0.651693,  -0.405118,
   -0.199855, -0.0102315,   0.180654,   0.389663,   0.642456,          1,          1,
           1,   0.712418,   0.526304,   0.388497,   0.278655,   0.186109,   0.103817,
   0.0263252, -0.0511668 >

数値を見ている限りは良く似ています。でもどれぐらい似ているのか分りにくいので、両方の差分をとってみましょう。

rslt[*,16] - rslt[*,14]
<  0.000675234,            0, -0.000638096,  -0.00116506,  -0.00152078,  -0.00166305,
   -0.00157029,  -0.00124316, -0.000705475,            0,            0,            0,
  3.81456e-005,  0.000109226,  0.000204127,  0.000291396,  0.000335864,  0.000308367,
   0.000193589,            0,            0,            0,  0.000565533,  0.000919033,
    0.00105694,  0.000996994,   0.00077539,  0.000440817, 4.89515e-005, -0.000342914 >

最大の差が -.0.00166305 です。この数字からすると上下のずれは 1% 位のようです。より詳しく上下の差の比率を計算して見ましょう。

vct=rslt[*,16], vct[1]=1, <<30,@|(rslt[_n,16] - rslt[_n,14])/vct[_n]>>
<      0.0114214,             0,     0.0106513,    0.00934579,    0.00766412,    0.00581638,
      0.00399282,    0.00234534,   0.000986351,             0,             0,             0,
   -5.85331e-005,  -0.000269614,   -0.00102137,    -0.0284804,    0.00185916,   0.000791368,
     0.000301327,             0,             0,             0,   0.000793822,     0.0017462,
      0.00272059,    0.00357787,    0.00416633,    0.00424608,    0.00185949,    0.00670188 >

rslt[*,16][1] == 0 で単純に割り算するとエラーになるので、ここだけ 1 に修正して差分の比率を計算しています。やはり 1% の誤差と見れるようです。

sf の便利さを分ってもらえますでしょう。ぜひとも御自分の手で色々と計算してみてください。

△φ = 0 計算結果の考察

△φ=0 は二階の偏微分方程式です。初期条件として一つの閉じた境界における初期値だけでは不十分です。その微分値も必要になります。

上のプログラムでは微分初期値を与えていません。暗黙のうちに無限遠点の値が 0 であるとの、二つ目の閉じた境界値が与えて解いています。でも無限遠点との二つの境界の値が与えられれば、全世界の値まで定まってしまうのは意外な気もします。三つ以上の閉じた境界値が与えられたら、それに合わせた △φ=0 の解が一意的に定まるのですから。

△φ=0 に浮いた導体を置いて見ます

(19,17),(19,18),(20,17),(20.18) の閉じた四点に浮いた金属を置いてみます。この四点が同じ電位になる条件を追加して △φ=0 偏微分方程式を解いてみます。

//@@
# coding=shift_jis
# 07.10.28 偏微分方程式 ラプラスの式 ≡ ∂^2φ/∂y^2 + ∂^2φ/∂x^2 == 0 電場
# 30 x 30  matrix position [15,10] と [15,20] の位置に 3 position 正方形を置き
# 0`V と 1`V の電位とする
# 外周 [:,0], [0,:], [:,N-1], [N-1,:] の位置は free:None とする
# <== free の位置の値は gradient より interpolate する
import sf
N = 30
dctBoundary = {}    # 境界条件を設定する辞書

# とりあえず 30x30 行列の全部を 計算可能:True にする
for index in sf.mrng(N,N):
    dctBoundary[index] = True

# x,y position [10, 15] にある 0`V 電位の設定
for index in sf.masq( (9,3), (14,3) ):
    dctBoundary[index] = -1

# x,y position [20,15] にある 1`V 電位の設定
for index in sf.masq( (19,3), (14,3) ):
    dctBoundary[index] = 1

# Gnd が取れていない導体をおく。境界は全て同じ電位になる
dctBoundary[17,19] = 'i'
dctBoundary[18,19] = 'i'
dctBoundary[17,20] = 'i'
dctBoundary[18,20] = 'i'

iterationAt = 40
for dctRslt, index in sf.itrSlvPDE(dctBoundary, iterationAt):
    j,k=index
    val = dctBoundary[index]
    if val == 'i':
        dctRslt[index] = dctRslt[16,19]
    else:
        dctRslt[j,k] = 1.0/4 *( dctRslt[j-1,k]
                              + dctRslt[j+1,k]
                              + dctRslt[j,k-1]
                              + dctRslt[j,k+1])

modelAt = sf.make3dRowCol(dctRslt)
sf.putSf(dctRslt, "rslt")
//@@@

dctBoundary[17,19] = [(16,19)] 以下の四行が浮いた導体の設定です。このメッシュ位置でで初期条件辞書の値を 'i' にしています。このメッシュ位置は dctRslt[16,19] の値と同じにしています。浮いた導体であるため、電位を [16,19] のメッシュ位置と同じにしたわけです。下のような結果になります。


もっともらしい 3D 図形がえられました。(17,19) の辺りに歪が見られます。こんなラフな計算でも、解の全体としての特徴は出てくるものです。

□ψ(x,t) == 0 : x,t 一次元波動方程式

∂^2ψ/∂x^2 - ∂^2ψ/∂t^2 == 0 の一次元波動方程式を sf.solvePDE(..) 関数で解きます。

∂^2ψ/∂x^2 - ∂^2ψ/∂t == 0 の漸化式

∂^2ψ/∂x^2, ∂^2ψ/∂t の二階微分は次のように離散化できます

∂^2ψ/∂x^2[j,k] ≒ (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 ------------ (1)
∂^2ψ/∂t^2[j,k] ≒ (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2 ------------ (2)
0 == ∂^2ψ/∂x^2 - ∂^2ψ/∂t == (1)式 - (2)式 の関係は次のように変形できます。
0 == (1)式 - (2)式 
  == (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 - (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2
∴
ψ[j,k+1] = (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])Δt^2/Δx^2 - ψ[j,k-1] + 2*ψ[j,k]
ψ[j,k] = 2*ψ[j,k-1] - ψ[j,k-2] + (ψ[j+1,k-1] + ψ[j-1,k-1] - 2*ψ[j,k-1])Δt^2/Δx^2 -------- (3)

ようは ψ[j,k+1] と時間を進めた値を、時間を進める前の周辺の値から導かせるということです。これは漸化式によって波形を掃きだしながら滑らかにしていく操作とも合致しています。

ここでは Δt^2 == Δx^2 を使った単純化をしていません。右辺の ψ[j,k] が打ち消しあってしまって嫌な感じがするためです。私の直感で Δt, Δx の区別を残しました。もしかしたら、単純化しても誤差は同じかもしれません。よく考えていません。Δt, Δx を残しても誤差が増えることにはなりません。プログラム・コードもさして変わりません。これでいきます。

∂^2ψ/∂x^2 - ∂^2ψ/∂t == 0 数値解

初期条件は左端で一周期だけ強制的に sin 波形で振動させます。その残りの時間は固定端とします。右端は自由端とします。時刻 0 における初期値は ψ(x,0) = 0 とします。ψ(x,1) = 0 の初期条件も (3) 式に k-2 の項があるので必要です。これは ∂ψ/∂t 初期値を与えることと等価です。

//@@
# coding=shift_jis
# 07.09.19 偏微分方程式 ∂^2ψ/∂x^2 - ∂^2ψ/∂t^2 == 0
#ψ[j,k] = (ψ[j+1,k-1] + ψ[j-1,k-1] - 2*ψ[j,k-1])Δt^2/Δx^2 - ψ[j,k-2] + 2*ψ[j,k-1] -------- (3)
# ただし Δt = 1/N, Δx^2 = 1/M
import sf
import math
N = 30      # t
timeWidth = 1.5
#timeWidth = 2.0
#N = 40      # t
M = 20      # x
dctBoundary = {}
for index in sf.mrng(M,N):
    dctBoundary[index] = True

# 時刻 0 における初期位置 ψ(x,0) = 0, ψ(x,1) = 0 を設定する
for i in range(M):
    dctBoundary[i,0] = 0
    dctBoundary[i,1] = 0

# 左端の時間軸の固定値を設定する。
for i in range(10):     # 10 まで sin 波駆動する
    dctBoundary[0,i] = math.sin(2*math.pi* i/10.0)

for i in range(10,N):   # 10 以降は 0 に固定
    dctBoundary[0,i] = 0

def fncPDE(dctRslt, index):
    j,k=index
    dctRslt[j,k] = 2*dctRslt[j,k-1] - dctRslt[j,k-2]\
                 + (timeWidth *M/N)**2*(   dctRslt[j+1,k-1]
                           + dctRslt[j-1,k-1]
                           - 2*dctRslt[j,k-1]
                         )

#import pprint as pp; print pp.pprint(dctBoundary)
dctRslt = sf.solvePDE(dctBoundary, fncPDE, 40)
modelAt = sf.make3dRowCol(dctRslt)

sf.putSf(dctRslt)
//@@@

プログラムを見直してみると漸化式が少し複雑になっています。あらためて偏微分方程式を解くことは、初期条件と漸化式を定めることだと思わされます。

計算結果は こんな具合です。



□ψ = 0 計算結果の考察

Δx == Δt にするのを嫌って少し複雑な漸化式にしましたが、Δx == Δt にしても、とくに計算結果は変わりませんでした。dctRslt[j,k-1] を打ち消し合わせた単純な漸化式を使っても良いようです。

自由端で倍近くマイナス方向に触れているのが引っかかりますが、それ以外は波動方程式として予測される振る舞いをしています。右側の自由端は線形外挿補間をしているのですが、それが右端でのマイナス方向へのオーバーシュートを発生させています。これは右端での微分値を与えているのだと解釈しています。

∂^2ψ/∂t^2 - ∂^2ψ/∂x^2 - ∂^2ψ/∂y^2 == 0 : 二次元波動方程式

∂^2ψ/∂t^2 - ∂^2ψ/∂x^2 - ∂^2ψ/∂y^2 == 0 の二次元波動方程式を sf.solvePDE(..) 関数で解きます。

漸化式は一次元の波動方程式を二次元に簡単に拡張できます。でも、その数値解は簡単に発散するようになります。時間変化に対し空間変化が小さくなるように刻み幅を注意深く選択することが必要になります。

∂ψ^2/∂t^2 - ∂ψ^2/∂x^2 - ∂ψ^2/∂y^2 == 0 の漸化式

∂^2ψ/∂x^2, ∂^2ψ/∂t の二階微分は次のように離散化できます
∂^2ψ/∂x^2[j,k,m] ≒ (ψ[j+1,k,m] + ψ[j-1,k,m] - 2*ψ[j,k,m])/Δx^2 ------------ (1)
∂^2ψ/∂y^2[j,k,m] ≒ (ψ[j,k+1,m] + ψ[j,k-1,m] - 2*ψ[j,k,m])/Δy^2 ------------ (2)
∂^2ψ/∂t^2[j,k,m] ≒ (ψ[j,k,m+1] + ψ[j,k,m-1] - 2*ψ[j,k,m])/Δt^2 ------------ (3)
0 == ∂^2ψ/∂x^2 + ∂^2ψ/∂y^2 - ∂^2ψ/∂t == (1)式 + (2)式 - (3)式 の関係は次のように変形できます。
0 == (1)式 + (2)式 - (3)式 
  == (ψ[j+1,k,m] + ψ[j-1,k,m] - 2*ψ[j,k,m])/Δx^2 
   + (ψ[j,k+1,m] + ψ[j,k-1,m] - 2*ψ[j,k,m])/Δy^2 
   - (ψ[j,k,m+1] + ψ[j,k,m-1] - 2*ψ[j,k,m])/Δt^2
∴
ψ[j,k,m+1] == 2*ψ[j,k,m] - ψ[j,k,m-1] 
             + (ψ[j+1,k,m] + ψ[j-1,k,m] - 2*ψ[j,k,m])Δt^2/Δx^2 
             + (ψ[j,k+1,m] + ψ[j,k-1,m] - 2*ψ[j,k,m])Δt^2/Δy^2 
             
∴
ψ[j,k,m] == 2*ψ[j,k,m-1] - ψ[j,k,m-2] 
             + (ψ[j+1,k,m-1] + ψ[j-1,k,m-1] - 2*ψ[j,k,m-1])Δt^2/Δx^2 
             + (ψ[j,k+1,m-1] + ψ[j,k-1,m-1] - 2*ψ[j,k,m-1])Δt^2/Δy^2 

最後の式の最初の二項、すなわち 2*ψ[j,k,m-1] - ψ[j,k,m-2] の項は時間に関して線形外挿補間をしている項と解釈できます。その後の項は周囲の空間部分からの影響を表していると理解できます。

二次元以上の波動方程式では、この周囲の空間部分からの影響が小さくなるような Δx, Δy, Δt を選択してやらないと簡単に数値解が発散します。一次元の場合は、そのような注意をしてやらなくても発散しませんでした。時間の変化にともなって空間の形を保ったままに平行移動していくだけだったからです。二次元以上では波形の形を変えながら移動していくため、漸化式を時間変化に対し空間変化を小さくするものにしてやらねばなりません。空間変化が大きくなると、それが漸化式の繰り返しにより増幅されていく状態に陥ります。そうなると数値解が発散します。

なお、漸化式のインデックスの増加させる順序も計算結果に大きく利いてきます。上の漸化式の時間 m で表しているインデックスは、下の python コードでは最初の要素にしてあります。空間のインデックス j, k で形を整えた後に時間軸での形を整える順序にしてやらないと計算結果が発散します。

∂^2ψ/∂t^2 - ∂^2ψ/∂x^2 - ∂^2ψ/∂x^2 == 0 数値解

初期条件は中央で 1 の 凸 領域としています。周囲は自由端にしています。sf.py 側の線形外挿補間にまかせています。

∂^2ψ/∂t^ の二回微分なので初期値は dctBoundary[0,?,?] の他に dctBoundary[1,?,?] も設定しています。

//@@
# 07.11.01 OK
import sf
import math
#N = 100
T = 200         # time mesh
N = 30          # x, y mesh
dt =0.1         # Δt
#dt =0.2         # Δt
#dt =1          #こっちを大きくすると発散する

dd =2.0         # Δx, Δy
#dd =0.2        # こっちを小さくすると発散する

dctBoundary = {}
for index in sf.mrng(T,N,N):
    dctBoundary[index] = True

# 時刻 0 における初期位置 ψ(x,y,0) = 0 を設定する
for index in sf.mrng(N,N):
    dctBoundary[(0,)+index]=0

# 中央が盛り上がった初期波形値とする
for index in sf.mrng([N/2-5,N/2+5],[N/2-5,N/2+5]):
    dctBoundary[(0,)+index] = 1

# 波動方程式は二階微分なので、時刻 1 における初期位置 ψ(x,y,1) = 0 も設定する
for j,k in sf.mrng(N,N):
    dctBoundary[1,j,k] = dctBoundary[0,j,k]

lstDbg = [2,2,2]
inCount = 0
for dctRslt, index in sf.itrSlvPDE(dctBoundary, 30):
    if tuple(lstDbg) == index:
        #import pdb; pdb.set_trace()
        print "loopCount:",inCount
        inCount += 1
    m, j,k = index
    dctRslt[m,j,k] = 2*dctRslt[m-1,j,k]-dctRslt[m-2,j,k]+dt*dt/(dd*dd)*(
                            dctRslt[m-1,j+1,k]+dctRslt[m-1,j-1,k]
                          + dctRslt[m-1,j,k+1]+dctRslt[m-2,j,k-1]
                          - 4*dctRslt[m-1,j,k] )
sf.putSf(dctRslt,"rslt")
//@@@
//copy \#####.### temp.py /y

計算結果は下のようになります。sf 式を使って三階の rslt.val テンソルから時間軸断面行列を抜き出して sfGrapy.py で 3D 表示させます

~sfGraph.py(rslt[10,*,*])
index:(14, 14)  max:     1
index:(11, 0)  min:-2.58675e-020

~sfGraph.py(rslt[60,*,*])
index:(14, 14)  max: 0.988547
index:(11, 11)  min:-0.157703

~sfGraph.py(rslt[120,*,*])
index:(14, 23)  max: 0.576194
index:(14, 14)  min:-0.504513

(□ - 1)ψ == 0 :Klein Goldon 方程式

□ψ == 0 を少しだけ変形した Klein Goldon 方程式 (∂^2/∂x^2 - ∂^2/∂t^2 - 1)ψ== 0 を数値解析的に解いてみましょう。

(∂^2/∂x^2 - ∂^2/∂t^2 -1)ψ== 0 の漸化式

∂^2ψ/∂x^2, ∂^2ψ/∂t の二階微分は次のように離散化できます

∂^2ψ/∂x^2[j,k] ≒ (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 ------------ (1)
∂^2ψ/∂t^2[j,k] ≒ (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2 ------------ (2)
0 == ∂^2ψ/∂t^2 - ∂^2ψ/∂x - ψ == (1)式 - (2)式 -ψ の関係は次のように変形できます。
0 == (1)式 - (2)式 - ψ
  == (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 - (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2 - ψ(j,k)
∴
ψ[j,k+1] = (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])Δt^2/Δx^2 - ψ[j,k-1] + 2*ψ[j,k] - ψ(j,k)Δt^2
ψ[j,k] = 2*ψ[j,k-1] - ψ[j,k-2] 
       + (ψ[j+1,k-1] + ψ[j-1,k-1] - 2*ψ[j,k-1])Δt^2/Δx^2 
       - ψ(j,k-1)Δt^2 -------- (3)

∂^2ψ/∂t^2 - ∂^2ψ/∂x -1 == 0 の数値解

時刻 0 において x 軸の中点に質点があるものとします。速度 は 0 とします。

//@@
# coding=shift_jis
# 07.10.18
# ψ[j,k] = 2*ψ[j,k-1] - ψ[j,k-2] 
#       + (ψ[j+1,k-1] + ψ[j-1,k-1] - 2*ψ[j,k-1])Δt^2/Δx^2 
#       - ψ(j,k-1)Δt^2

import sf

timeWidth = 0.2
N = 30      # t
#N = 40      # t
M = 20      # x
dctBoundary = {}

for index in sf.mrng(M,N):
    dctBoundary[index] = True

#                                   ┌──┐
# 時刻 0 における初期位置 ψ(0,x) ─┘    └─ を設定する
for i in range(M):
    if (i == int(M/2) ):
        dctBoundary[i,0] = 1
        dctBoundary[i,1] = 1    # 二階微分なので、初期条件を もう一つ追加する
    else:
        dctBoundary[i,0] = 0
        dctBoundary[i,1] = 0    # 二階微分なので、初期条件を もう一つ追加する

deltaT = timeWidth/N
deltaX = 1.0/M

def fncPDE(dctRslt, index):
    j, k = index
    dctRslt[j,k] = (  
            2*dctRslt[j,k-1]- dctRslt[j,k-2]
          + (dctRslt[j+1,k-1]- dctRslt[j-1,k-1]- 2*dctRslt[j,k-1])*deltaT**2/deltaX**2
          - dctRslt[j,k-1]*deltaT**2
    )

dctRslt = sf.solvePDE(dctBoundary, fncPDE, 30)
modelAt = sf.make3dRowCol(dctRslt)
sf.putSf(dctRslt,"rslt")
//@@@

計算結果は こんな具合です。物理的な直感が予想するように、振動しながら広がっていく粒子になっています。時間の進行に伴って発散しています。timeWidth を大きくすると発散が顕著になります。この理由が計算誤差の累積なのか Klein Goldon 方程式の本質なのか、漸化式の導出に誤りがあるのか、私には解かりません。(deltaT を変えても発散のしかたは変わりません。)御存知の方はお教えください。





∂^2ψ/∂x^2 - ∂^2ψ/∂t -1 == 0 :Klein Goldon もどき方程式

□ψ - 1 == 0 を数値解析的に解いてみましょう。実は、こっちを Klein Goldon 方程式だと最初は勘違いしていました。本当の Klein Goldon 方程式との比較の意味で価値があるかと考え、残します。

∂^2ψ/∂x^2 - ∂^2ψ/∂t -1 == 0 の漸化式

∂^2ψ/∂x^2, ∂^2ψ/∂t の二階微分は次のように離散化できます

∂^2ψ/∂x^2[j,k] ≒ (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 ------------ (1)
∂^2ψ/∂t^2[j,k] ≒ (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2 ------------ (2)
0 == ∂^2ψ/∂t^2 - ∂^2ψ/∂x - 1== (1)式 - (2)式 -1 の関係は次のように変形できます。
0 == (1)式 - (2)式 - 1
  == (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])/Δx^2 - (ψ[j,k+1] + ψ[j,k-1] - 2*ψ[j,k])/Δt^2 - 1
∴
ψ[j,k+1] = (ψ[j+1,k] + ψ[j-1,k] - 2*ψ[j,k])Δt^2/Δx^2 - Δt^2 - ψ[j,k-1] + 2*ψ[j,k]
ψ[j,k] = 2*ψ[j,k-1] - ψ[j,k-2] -Δt^2 + (ψ[j+1,k-1] + ψ[j-1,k-1] - 2*ψ[j,k-1])Δt^2/Δx^2 -------- (3)

∂^2ψ/∂t^2 - ∂^2ψ/∂x -1 == 0 の数値解

時刻 0 において x 軸の中点に質点があるものとします。その速度 は 0 とします。

//@@
# coding=shift_jis
import sf

timeWidth = 0.2
#N = 30      # t
N = 40      # t
M = 20      # x
dctBoundary = {}

for index in sf.mrng(M,N):
    dctBoundary[index] = True

#                                   ┌──┐
# 時刻 0 における初期位置 ψ(0,x) ─┘    └─ を設定する
for i in range(M):
    if (i == int(M/2) ):
        dctBoundary[i,0] = 1
        dctBoundary[i,1] = 1    # 二階微分なので、初期条件を もう一つ追加する
    else:
        dctBoundary[i,0] = 0
        dctBoundary[i,1] = 0    # 二階微分なので、初期条件を もう一つ追加する

def fncPDE(dctRslt, index):
    if index == (10,1):  # to debug
        print index
    j, k = index
    dctRslt[j,k] = 2*dctRslt[j,k-1] - dctRslt[j,k-2]\
                 - 1*(timeWidth/N)**2\
                 + (timeWidth *M/N)**2*(   dctRslt[j+1,k-1]
                           + dctRslt[j-1,k-1]
                           - 2*dctRslt[j,k-1]
                   )


dctRslt = sf.solvePDE(dctBoundary, fncPDE, 100)
modelAt = sf.make3dRowCol(dctRslt)
sf.putSf(dctRslt,"rslt")
//@@@

計算結果は こんな具合です。本物とはことなり左右対称になっています。



∂f(x,t)/∂x)^2 + ∂f(x,t)/∂t=g(q) Hamilton Jacobi 偏微分方程式 1

解析力学の教科書で、調和振動子の Hamilton Jacobi 偏微分方程式 (∂f(x,t)/∂x)^2 + ∂f(x,t)/∂t=g(q) の解は f(q,a) - a t の変数分離型になると書いてあります。これがずっと納得できませんでした。sf.solePDE(..) で Hamilton Jacobi 方程式の数値解を求めて、この変数分離の可能性を数値実験によって確認してみたいと思います。

( f(x,t)/∂x)^2 + ∂f(x,t)/∂t==0 の漸化式

単純な g(q) = 0 のときの ( f(x,t)/∂x)^2 + ∂f(x,t)/∂t==0 に話を限定します。まず、下の微分式を前提とします。

∂f/∂t ≒ (f(j,k) - f(j-1,k))/Δt ------------ (1) 式
∂f/∂x ≒ (f(j-1,k) - f(j-1,k-1))/Δx -------- (2) 式

0 == (2)式^2 + (1)式 の関係式から下のように漸化式を導きます。

0 == (2)式^2 + (1)式 
0 == [(f(j-1,k) - f(j-1,k-1))/Δx]^2 + (f(j,k) - f(j-1,k))/Δt
0 == [(f(j-1,k) - f(j-1,k-1))]^2 Δt/Δx^2 + (f(j,k) - f(j-1,k))
∴
f(j,k) == f(j-1,k) - [(f(j-1,k) - f(j-1,k-1))]^2 Δt/Δx^2 

( f(x,t)/∂x)^2 + ∂f(x,t)/∂t==0 の数値解

//@@
# coding=shift_jis
# 07.10.01 偏微分方程式 ( f(x,t)/∂x)^2 + x^2 +∂f(x,t)/∂t==0 OK
# 
# 初期条件 エネルギー一定条件 p^2 +x^2 == const     ×f(x,0) = sin(2πx)
# f(t,x):縦軸を時間軸、横軸を x 軸にする matrix[j,k]: j は時間、k は位置の x 座標
# f(j,k) == f(j-1,k) - [(f(j-1,k) - f(j-1,k-1))]^2 Δt/Δx^2 - x^2 Δt
# ただし Δt = 1/N, Δx^2 = 1/M
# 初期条件 f(0,x) = sin(2π x)
import sf
import math
#rate = 1.0
#timeSpan = 2.0  # timeSpan = 1 のままだと、t 軸メッシュの細かさを指定するだけになる
timeSpan = 1.3
#timeSpan = 1.0  # timeSpan = 1 のままだと、t 軸メッシュの細かさを指定するだけになる
#timeSpan = 0.5  # timeSpan = 1 のままだと、t 軸メッシュの細かさを指定するだけになる
N = 40      # t
#N = 20      # t
M = 20      # x
#M = 40      # x
dctBoundary = {}
for index in sf.mrng(M,N):
    dctBoundary[index] = True

# 初期条件 f(x,0) = 1
for i in range(M):     # 10 まで駆動する
    #dctBoundary[i,0] = 0.013*math.sin( 2*math.pi*i/(M-1) ) # 発散を始めます
    #dctBoundary[i,0] = 0.010*math.sin( 2*math.pi*i/(M-1) ) # 発散を始めます
    dctBoundary[i,0] = 0.008*math.sin( 2*math.pi*i/(M-1) ) # 発散を始めます
    #dctBoundary[i,0] = 0.002*math.sin( 2*math.pi*i/(M-1) ) # W(q,a)≡Wl(q)Wr(a) に限定するときは a が虚数になる
    #dctBoundary[i,0] = (2.0 - (1.0*i/M)**2 )**0.5
    #dctBoundary[i,0] = 0.1*i/M
    #dctBoundary[i,0] = 0.1*math.asin(1.0*i/M) + 0.01*1.0*i/M*(1-(1.0*i/M)**2)**0.5
    #dctBoundary[i,0] = 0.5*math.asin(1.0*i/M) + 0.5*1.0*i/M*(1-(1.0*i/M)**2)**0.5 # too large

# dctRslt を繰り返し計算し収束を待つ。収束はグラフを目視することで判定する
deltaX = 1.0/M
deltaT = timeSpan/N
def fncPDE(dctRslt, index):
    j,k = index
    dctRslt[j,k] = dctRslt[j,k-1]\
                  - ( dctRslt[j,k-1] - dctRslt[j-1,k-1] )**2 *deltaT/deltaX**2\
                  - (deltaX*k/M)**2*deltaT


dctRslt = sf.solvePDE(dctBoundary, fncPDE, inIteration = 40)
modelAt = sf.make3dRowCol(dctRslt)
sf.putSf(dctRslt, "rslt")
//@@@

計算結果は こんな具合です。


この漸化式は直ぐに発散する性質があります。dctRslt[j,k-1] - dctRslt[j-1,k-1] )**2 *deltaT/deltaX**2 の項があり、微分値が大きくなる、即ち変移が激しくなると、正帰還がかかり発散してしまいます。解ける初期条件が限られます。

それでも、発散しない範囲で色々な初期条件を与えてみると、偏微分方程式の解が変数分離型にならないものもあることが分ります。解析力学の Hamilton Jacobi の説明に書いてあった説明は、「変数分離型になる」のではなく「変数分離型にできる解を選ぶ」の意味であることが分りました。

しかし まあ sf.solvePDE(..) は、よくもこんな偏微分方程式まで解いてくれるものです。我ながら感心します。

div E = 0 : ベクトル偏微分方程式

こんどはベクトル値の偏微分方程式 div E を解いてみましょう。

div vctE == 0 の漸化式

div vctE = 0:ベクトル偏微分方程式を離散化してみましょう。下で vctE はベクトル値です。スカラー値ではありません。
        ∂vctE/∂x ≒ (vctE(j,k) - vctE(j-1,k))/Δx ------------ (1) 式
        ∂vctE/∂y ≒ (vctE(j,k) - vctE(j,k-1))/Δy ------------ (2) 式
        ∂vctE/∂x ≒ (vctE(j+1,k) - vctE(j,k))/Δx ------------ (3) 式
        ∂vctE/∂y ≒ (vctE(j,k+1) - vctE(j,k))/Δy ------------ (4) 式

これを組み合わせて div vctE に対して、下のような漸化式を作れます。(注意:この漸化式は十分条件の漸化式にすぎません。元々 div vctE = 0 だけでは一意的に解は定まりません。でも この漸化式で div vctE == 0 を満たす特別な解を作ることはできます。)

0 == [ (1) 式 + (2) 式 ] - [ (3) 式 + (4) 式 ]
  == (vctE(j,k) - vctE(j-1,k))/Δx
     (vctE(j,k) - vctE(j,k-1))/Δy
    -(vctE(j+1,k) - vctE(j,k))/Δx
    -(vctE(j,k+1) - vctE(j,k))/Δy
  == 2vctE(j,k)(1/Δx + 1/Δy) - [ vctE(j-1,k)+vctE(j+1,k)]/Δx - [vctE(j,k-1)+vctE(j,k+1)]/Δy
  == vctE(j,k)2(Δx + Δy)/(ΔxΔy) - [ vctE(j-1,k)+vctE(j+1,k)]/Δx - [vctE(j,k-1)+vctE(j,k+1)]/Δy
∴
  vctE(j,k)= 0.5/(Δx + Δy){ [vctE(j-1,k)+vctE(j+1,k)]Δy + [vctE(j,k-1)+vctE(j,k+1)]Δx }

下のような漸化式を作ることもできます。でも この漸化式は左下から上側に吐き出していくだけであり、全体の形を修練させることができないので、上の四方の周囲を平均化する漸化式を使います。

∂vctE/∂x ≒ (vctE(j,k) - vctE(j-1,k))/Δx ------------ (1) 式
∂vctE/∂y ≒ (vctE(j,k) - vctE(j,k-1))/Δy ------------ (2) 式

0 == (1) 式 + (2) 式
  == (vctE(j,k) - vctE(j-1,k))/Δx + (vctE(j,k) - vctE(j,k-1))/Δy
  == vctE(j,k)(1/Δx + 1/Δy)- vctE(j-1,k))/Δx - vctE(j,k-1))/Δy
∴
vctE(j,k) = 1/(Δx +Δy)[ vctE(j-1,k))Δy + vctE(j,k-1)Δx]

div vctE == 0 の数値解

ベクトルの漸化式であっても sf.solvePDE(.) は処理できます。ただしデフォルト値を scalar の 0 にできないので dfltValue=vs.array([0,0]) のようにデフォルト・ベクトル値を明示的にユーザーが設定する必要があります。

//@@
# coding=shift_jis
# 07.10.08  
# test div E = 0    NG
import sf
import visual as vs

N = 30
M = 30
dctBoundary = {}    # 境界条件を設定する辞書

# とりあえず 30x30 行列の全部を 計算可能:True にする

for index in sf.mrng(N,M):
    dctBoundary[index] = True

scale = 0.1
r=5.0*scale
# matrix position [15,15] を中心とする半径 5 の球表面の電場ベクトルの設定
#N=100
#lstAt =[(r,0)]
#for theta in sf.arSqnc(0,N, 2*math.pi/N):
#    y,x = math.floor(r*math.cos(theta)), math.floor(r*math.sin(theta))
#    if lstAt[-1] != (y,x):
#        lstAt.append((y,x))
#print lstAt
#             (15,20)              (x,y) 座標 
#       (14,19)     (16,19)
#    (12,18)           (17,18)
#  (11,17)                (19,17)
# (10,16)                  (20,16)
# (10,15)                  (20,15)
# (10,14)                  (20,14)
#  (11,13)                (18,13)
#    (12,12),           (17,12)
#       (14,11),    (16,11)
#             (15,10)

posLeftAt=vs.array([15,20])
posCenterAt = vs.array([15,15])
dctBoundary[tuple(posLeftAt)] = (posLeftAt - posCenterAt)*scale/r**3

posLeftAt=vs.array([14,19])
posRightAt=vs.array([16,19])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([12,18])
posRightAt=vs.array([17,18])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([11,17])
posRightAt=vs.array([19,17])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,16])
posRightAt=vs.array([20,16])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,15])
posRightAt=vs.array([20,15])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,14])
posRightAt=vs.array([20,14])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([11,13])
posRightAt=vs.array([18,13])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([12,12])
posRightAt=vs.array([17,12])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([14,11])
posRightAt=vs.array([16,11])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([15,10])
dctBoundary[tuple(posLeftAt)] = (posLeftAt - posCenterAt)*scale/r**3

deltaY = scale
deltaX = scale
def fncPDE(dctRslt, index):
    j,k=index
    dctRslt[j,k]= 0.5/(deltaX + deltaY)*(
                            (dctRslt[j-1,k]+dctRslt[j+1,k] )*deltaY
                          + (dctRslt[j,k-1]+dctRslt[j,k+1] )*deltaX
                  )

dctRslt = sf.solvePDE(dctBoundary, fncPDE, inIteration=40, dfltValue=vs.array([0,0]))
sf.putSf(dctRslt, "rslt")

import visual as vs
for index in sf.mrng(N,N):
    vs.arrow(pos=vs.vector(index)-[15,15], axis =dctRslt[index]
            , shaftwidth = 0.3, color = vs.color.red)
//@@@

円形の初期値を設定したので、dctBoundary の設定に手間がかかっています。でも sf.sovePDE(..) の使い方は dfltValue 引数が増えた以外は、今までと同じです。漸化式の中身がベクトルの加減乗除算になっているのですが、Duck Typing の python コードでは、スカラーのときと同じに見えてしまいます。同じように扱えてしまいます。

計算結果の可視化については、今までのように sf.make3dRowCol(..) 関数に任せるわけにはいきません。なんせベクトル値の分布なのですから。

ここでは vpython の arrow(..) オブジェクトを使って dctRslt[,] 行列の全ての要素をベクトルを表示させています。下の結果になりました。


div E == 0 の方程式では閉じた境界での表面のベクトル分布を定めれば、空間全体のベクトル分布が定まるとなった。驚きだ....ちっょとまて。そんな馬鹿な! 電磁気の Lorentz 条件式:∂φ/∂t + ∂A/∂t == 0 が成り立つとき、閉じた境界表面での <φ,A >分布を定めるだけで空間全体の <φ,A> 分布が定まるか? 連続方程式 :∂ρ/∂t + ∂j/∂t == 0 が成り立つとき、閉じた境界表面での <ρ,j>分布を定めるだけで空間全体の <ρ,A> 分布が定まるか?そんな訳はない。

そもそも ∂Ex/∂x + ∂Ey/∂y == 0 の式には、<Ex,Ey> と二成分の自由度がある。それなのに div E == 0 の一つの式だけで解が一意的に定まるとするのに無理がある。rot E:∂Ex/∂x - ∂Ey/∂y == 0 のような、もう一つの条件式を定めないと <Ex,Ey> 分布は決まらないだろう。

それでも上で作った漸化式は間違いではない。div E == 0 を満たす解であることには間違いない。ただし十分条件に過ぎなかったのだ。上で作り上げた漸化式の必要十分条件となるもう一つのベクトル偏微分方程式があるはずだ。

div vctE == 0, rot vctE == 0 の漸化式

電場の偏微分方程式ならば div vctE == 0 だけでなく rot vctE == 0 も満たすはずです。二つ目の偏微分方程式を与えれば、二成分の vctE も一意に定まるはずです。この二つを満たす漸化式を下のように求められます。

∂vctE[0]/∂x + ∂vctE[1]/∂y == 0:div
(vctE(j,k)[0] - vctE(j-1,k)[0]) + (vctE(j,k)[1] - vctE(j,k-1)[1]) == 0 --- (1)
 vctE(j,k)[0] + vctE(j,k)[1] == vctE(j-1,k)[0] + vctE(j,k-1)[1] ---------- (1a)
(vctE(j+1,k)[0] - vctE(j,k)[0]) + (vctE(j,k+1)[1] - vctE(j,k)[1]) == 0 --- (1')
 vctE(j,k)[0] + vctE(j,k)[1] ==  vctE(j+1,k)[0] + vctE(j,k+1)[1] --------- (1a')


∂vctE[1]/∂x - ∂vctE[0]/∂y == 0:rot
(vctE(j,k)[1] - vctE(j-1,k)[1]) - (vctE(j,k)[0] - vctE(j,k-1)[0]) == 0 --- (2)
-vctE(j,k)[0] + vctE(j,k)[1]  ==-vctE(j,k-1)[0] + vctE(j-1,k)[1] --------- (2a)

(vctE(j+1,k)[1] - vctE(j,k)[1]) - (vctE(j,k+1)[0] - vctE(j,k)[0]) == 0 --- (2')
 vctE(j,k)[0] - vctE(j,k)[1]  == vctE(j,k+1)[0] - vctE(j+1,k)[1] --------- (2a')

(1a)+(2a')
 vctE(j,k)[0] + vctE(j,k)[1] == vctE(j-1,k)[0] + vctE(j,k-1)[1] ---------- (1a)
 vctE(j,k)[0] - vctE(j,k)[1]  == vctE(j,k+1)[0] - vctE(j+1,k)[1] --------- (2a')
2vctE(j,k)[0] == vctE(j-1,k)[0] + vctE(j,k-1)[1] + vctE(j,k+1)[0] - vctE(j+1,k)[1] -- (1a)+(2a')

(1a')+(2a)
 vctE(j,k)[0] + vctE(j,k)[1] ==  vctE(j+1,k)[0] + vctE(j,k+1)[1] --------- (1a')
 vctE(j,k)[0] - vctE(j,k)[1]  == vctE(j,k+1)[0] - vctE(j+1,k)[1] --------- (2a')
2vctE(j,k)[1] == vctE(j+1,k)[0] + vctE(j,k+1)[1] - vctE(j,k-1)[0] + vctE(j-1,k)[1] -- (1a')+(2a)
<== 位置が symetric でない
(1a)-(2a')
 vctE(j,k)[0] + vctE(j,k)[1] == vctE(j-1,k)[0] + vctE(j,k-1)[1] ---------- (1a)
 vctE(j,k)[0] - vctE(j,k)[1]  == vctE(j,k+1)[0] - vctE(j+1,k)[1] --------- (2a')
2vctE(j,k)[1] == vctE(j-1,k)[0] + vctE(j,k-1)[1] - vctE(j,k+1)[0] + vctE(j+1,k)[1] -- (1a)-(2a')

(1a')-(2a)
 vctE(j,k)[0] + vctE(j,k)[1] ==  vctE(j+1,k)[0] + vctE(j,k+1)[1] --------- (1a')
-vctE(j,k)[0] + vctE(j,k)[1]  ==-vctE(j,k-1)[0] + vctE(j-1,k)[1] --------- (2a)
2vctE(j,k)[0] == vctE(j+1,k)[0] + vctE(j,k+1)[1]+ vctE(j,k-1)[0] - vctE(j-1,k)[1]  -- (1a')-(2a)

(1a)+(2a')+(1a')-(2a)
2vctE(j,k)[0] == vctE(j-1,k)[0] + vctE(j,k-1)[1] + vctE(j,k+1)[0] - vctE(j+1,k)[1] -- (1a)+(2a')
2vctE(j,k)[0] == vctE(j+1,k)[0] + vctE(j,k+1)[1]+ vctE(j,k-1)[0] - vctE(j-1,k)[1]  -- (1a')-(2a)
4vctE(j,k)[0] == vctE(j-1,k)[0] + vctE(j,k-1)[1] + vctE(j,k+1)[0] - vctE(j+1,k)[1]
               + vctE(j+1,k)[0] + vctE(j,k+1)[1] + vctE(j,k-1)[0] - vctE(j-1,k)[1]

(1a')+(2a)+(1a)-(2a')
2vctE(j,k)[1] == vctE(j+1,k)[0] + vctE(j,k+1)[1] - vctE(j,k-1)[0] + vctE(j-1,k)[1] -- (1a')+(2a)
2vctE(j,k)[1] == vctE(j-1,k)[0] + vctE(j,k-1)[1] - vctE(j,k+1)[0] + vctE(j+1,k)[1] -- (1a)-(2a')
4vctE(j,k)[1] == vctE(j+1,k)[0] + vctE(j,k+1)[1] - vctE(j,k-1)[0] + vctE(j-1,k)[1]
               + vctE(j-1,k)[0] + vctE(j,k-1)[1] - vctE(j,k+1)[0] + vctE(j+1,k)[1]

(1a)+(2a'), (1a')+(2a) も div vctE == 0 だけでなく rot vctE == 0 を満たすのですが、式が対称でありません。これを使うと訳の解からない結果になります。対称性のより良い (1a)+(2a')+(1a')-(2a), (1a')+(2a)+(1a)-(2a') 式を使います。

div vctE == 0, rot vctE == 0 の数値解

初期条件は前回と一緒です。漸化式のみを入れ替えて数値解を求めてみます。

//@@
# coding=shift_jis
# 07.10.14  
# test div E = 0    NG
import sf
import visual as vs

N = 30
M = 30
dctBoundary = {}    # 境界条件を設定する辞書

# とりあえず 30x30 行列の全部を 計算可能:True にする

for index in sf.mrng(N,M):
    dctBoundary[index] = True

scale = 0.1
r=5.0*scale
# matrix position [15,15] を中心とする半径 5 の球表面の電場ベクトルの設定
#N=100
#lstAt =[(r,0)]
#for theta in sf.arSqnc(0,N, 2*math.pi/N):
#    y,x = math.floor(r*math.cos(theta)), math.floor(r*math.sin(theta))
#    if lstAt[-1] != (y,x):
#        lstAt.append((y,x))
#print lstAt
#             (15,20)              (x,y) 座標 
#       (14,19)     (16,19)
#    (12,18)           (17,18)
#  (11,17)                (19,17)
# (10,16)                  (20,16)
# (10,15)                  (20,15)
# (10,14)                  (20,14)
#  (11,13)                (18,13)
#    (12,12),           (17,12)
#       (14,11),    (16,11)
#             (15,10)

posLeftAt=vs.array([15,20])
posCenterAt = vs.array([15,15])
dctBoundary[tuple(posLeftAt)] = (posLeftAt - posCenterAt)*scale/r**3

posLeftAt=vs.array([14,19])
posRightAt=vs.array([16,19])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([12,18])
posRightAt=vs.array([17,18])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([11,17])
posRightAt=vs.array([19,17])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,16])
posRightAt=vs.array([20,16])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,15])
posRightAt=vs.array([20,15])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([10,14])
posRightAt=vs.array([20,14])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([11,13])
posRightAt=vs.array([18,13])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([12,12])
posRightAt=vs.array([17,12])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([14,11])
posRightAt=vs.array([16,11])
for x in range(posLeftAt[0], posRightAt[0]+1):
    posAt = vs.array([ x,posLeftAt[1] ])
    if (x == posLeftAt[0]) or (x == posRightAt[0]):
        dctBoundary[tuple(posAt)] = (posAt - posCenterAt)*scale/r**3
    else:
        dctBoundary[tuple(posAt)] = False

posLeftAt=vs.array([15,10])
dctBoundary[tuple(posLeftAt)] = (posLeftAt - posCenterAt)*scale/r**3

#lstAt =[ 1, 1]  #to debug
def fncPDE(dctRslt, index):
    j,k=index
    dctRslt[j,k][0]  = 0.25*( dctRslt[j-1,k][0] 
                            + dctRslt[j,k-1][1] 
                            + dctRslt[j,k+1][0] 
                            - dctRslt[j+1,k][1]
                            + dctRslt[j+1,k][0] 
                            + dctRslt[j,k+1][1] 
                            + dctRslt[j,k-1][0] 
                            - dctRslt[j-1,k][1]
                            )

    dctRslt[j,k][1]  = 0.25*( dctRslt[j+1,k][0] 
                            + dctRslt[j,k+1][1] 
                            - dctRslt[j,k-1][0] 
                            + dctRslt[j-1,k][1]
                            + dctRslt[j-1,k][0] 
                            + dctRslt[j,k-1][1] 
                            - dctRslt[j,k+1][0] 
                            + dctRslt[j+1,k][1]
                            )

dctRslt = sf.solvePDE(dctBoundary, fncPDE, inIteration=40, dfltValue=sf.sc.array([0.0, 0.0]) )
print dctRslt[10,17]
sf.putSf(dctRslt, "rslt")

import visual as vs
for index in sf.mrng(N,N):
    vs.arrow(pos=vs.vector(index)-[15,15], axis =dctRslt[index]
            , shaftwidth = 0.3, color = vs.color.red)
//@@@

数値解は下のようなりました。


微妙にバイアスがかかった解です。本来ならば電場分布は円対称にならねばなりません。漸化式の対称性の不十分さに引っ張られています。たぶん [j+1,k+1] などの項を追加した、より対称性の高い漸化式を使うべきでしょう。その漸化式が、先の div E で早とちりで計算したベクトル式の漸化式の場合に一致するのだと推測されます。でも未だ より高い対称性を持った漸化式を導けていません。どなたか御存知の方は御教えください。

視点を変えます。 rot E:∂Ex/∂x - ∂Ey/∂y == 0 条件を追加したら複素解析関数になってしまいます。ベクトルの方向が複素数にするとき違ってきますが。より対称性の高いの漸化式を作れたら、閉じていなくても、一部の曲線上での複素数値分布を定めるだけで複素空間全体での分布が定まることを示せるはずです。それは多価関数になること、即ちリーマン面の発生もありうるはずです。面白く遊べそうです。この意味でも良い漸化式を御存知の方は、それを御教えください。

(∂S/∂x)^2 + (∂S/∂y)^2 - K == 0 Hamilton Jacobi 偏微分方程式 2

Hamilton Jacobi 偏微分方程式 1 では変数分離解に限定して考えることをのべました。ですので、時間を含まない空間偏微分方程式を考えます。まだ空間に分布するポテンシャルまでは考えません。

先に結論を言うと、この偏微分方程式は性質の悪い偏微分方程式です。 △ψ = 0, □φ=0 のように素直な方程式ではありません。まず方程式の形が線形ではありません。初期条件により、解の形は大きく変わります。数値解を求めようとしても簡単に複素数値解に入り込んでしまいます。複素数値解に入り込ませないために許される初期条件は強く限定されます。

Hamilton Jacobi 方程式の場合は、偏微分方程式を解くというより、物理条件などから許される初期条件/境界条件を探し、それに合わせた漸化式を作ってやって解くという作業になってしまうと思います。闇雲に漸化式を作って数値計算するだけでは、複素数値解に入り込んだり発散したりするだけです。

(∂S/∂x)^2 + (∂S/∂y)^2 - K == 0 の初期条件と漸化式

S(x,y) == (K/2)^0.5x + K/2^0.5.y が、掲題の Hamilton Jacobi 方程式の一つの解です。(0,0) 位置の近傍で、この解方程式の値をの初期条件値を設定するものとし、この初期条件に合わせた漸化式を考えます。

まず下のように微分式を離散化します。

∂f/∂x ≒ (f(j,k) - f(j-1,k))/Δx ------------ (1) 式
∂f/∂y ≒ (f(j,k) - f(j,k-1))/Δy ------------ (2) 式
h == Δx == Δy とし Hamilton Jacobi 方程式 (1 式)^2 + (2 式)^2 - 1 == 0 より下のように漸化式を導けます。
0 == [f(j,k) - f(j-1,k)]^2 + [f(j,k) - f(j,k-1)]^2 + h^2 C
  == 2f(j,k)^2 [- 2f(j,k)f(j-1,k) + f(j-1,k)^2 + [- 2(f(j,k)f(j,k-1)+f(j,k-1)^2 ] + h^2 C
  == 2f(j,k)^2 - 2f(j,k)[f(j-1,k) + f(j,k-1)]+ [f(j-1,k)^2 + f(j,k-1)^2] + h^2 C

∴
0 ==  f(j,k)^2 -  f(j,k)[f(j-1,k) + f(j,k-1)]+ 0.5[f(j-1,k)^2 + f(j,k-1)^2 + h^2 K]
上の式を f(j,k) についての二次方程式とみなすと係数 b, c が次のように定まります。
b ≡ f(j-1,k) + f(j,k-1)
c ≡ 0.5[f(j-1,k)^2 + f(j,k-1)^2 + h^2 K]

この係数より f(j,k) が下のように定まります。
f(j,k) == 0.5*(-b + (b^2-4*c)^0.5)

上の漸化式は、x y 座標の原点 (0,0) の位置から右上に掃きだしていく形式の漸化式となっています。初期値の位置を (5,5) などに位置にとったら、(0,0) の値は定まらなくなります。

四方の周囲のの値に対して平均をとると考えて下のように漸化式を作ることもできます。でも、こちらの漸化式では複素数値解になるばかりでした。

0 == f(j,k)^2 - 0.5[f(j-1,k) + f(j,k-1) + f(j+1,k) + f(j,k+1)] f(j,k)
    + 0.25[f(j-1,k)^2 + f(j,k-1)^2 + f(j+1,k)^2 + f(j,k+1)^2] 
   - 0.5h^2 K

(∂S/∂x)^2 + (∂S/∂y)^2 - K == 0 の数値解

初期値は for pos in sf.masq([0+0j,3],[0+0j,3]):dctBoundary[pos] = sf.sc.array([pos[0]*scale/M +pos[1]*scale/N, pos[0]*scale/M +pos[1]*scale/N]) で与えています。

//@@
# coding=shift_jis
# 07.10.29
# test (∂S/∂x)^2 + (∂S/∂y)^2 - K = 0 Hamilton-Jacobi 2
# とりあえず K = 1
# 複素数値関数として扱う。二値の多価関数として扱う。そのため長さ 2 のベクトル値関数とする

import sf
import visual as vs

M = 30
#M = 50
N = M
dctBoundary = {}    # 境界条件を設定する辞書

# とりあえず 30x30 行列の全部を 計算可能:True にする

for index in sf.mrng(M,N):
    dctBoundary[index] = True

#scale = 1
#scale = 2
#scale = 0.2
scale = 0.02
# [0,0,][0,1][0,2][1,0][1,1][1,2][2,0[2,1][2,2] の位置で
# S(x,y) == (K/2)^0.5x + K/2^0.5.y の初期値だとします。 for pos in sf.masq([0,3],[0,3]): #dctBoundary[pos] = 1.0 #dctBoundary[pos] = 0.1 dctBoundary[pos] = sf.sc.array([pos[0]*scale/M +pos[1]*scale/N + 0j ,pos[0]*scale/M +pos[1]*scale/N + 0j]) #,-pos[0]*scale/M +pos[1]*scale/N]) #dctBoundary[pos] = 0.001 for n in range(3,N): dctBoundary[0,n] = 'i' for m in range(3,M): dctBoundary[m,0,] = 'j' lstAt =[ 1, 1, 3] #to debug h = float(scale)/M K = 1 #× K = -1.0 #K = 10 #K = 20 #K = 1000 import math def fncPDE(dctRslt, index): j,k=index if (dctBoundary[j,k] == 'i'): dctRslt[j,k] = 2*dctRslt[j,k-1] - dctRslt[j,k-2] return elif (dctBoundary[j,k] == 'j'): dctRslt[j,k] = 2*dctRslt[j-1,k] - dctRslt[j-2,k] return b = -( dctRslt[j-1,k][0] + dctRslt[j,k-1][0] ) c = 0.5*( dctRslt[j-1,k][0]**2 + dctRslt[j,k-1][0]**2 - h**2 * K ) cplxP = 0.5*(-b + (b**2-4*c)**0.5) cplxM = 0.5*(-b - (b**2-4*c)**0.5) if cplxP.imag < cplxM.imag: dctRslt[j,k][0] = cplxP else: dctRslt[j,k][0] = cplxM b = -( dctRslt[j-1,k][1] + dctRslt[j,k-1][1] ) c = 0.5*( dctRslt[j-1,k][1]**2 + dctRslt[j,k-1][1]**2 - h**2 * K ) cplxP = 0.5*(-b + (b**2-4*c)**0.5) cplxM = 0.5*(-b - (b**2-4*c)**0.5) if cplxP.imag >= cplxM.imag: dctRslt[j,k][1] = cplxP else: dctRslt[j,k][1] = cplxM dctRslt = sf.solvePDE(dctBoundary, fncPDE, inIteration=40, dfltValue=sf.sc.array([0+0j ,0+0j])) sf.putSf(dctRslt, "rslt") dctVw = {} for index in sf.mrng(M,N): #dctVw[index] = abs(dctRslt[index][0]) dctVw[index] = abs(dctRslt[index][1]) sf.make3dRowCol(dctVw) //@@@

数値解は下のようなりました。微妙に湾曲していますが、収束しきれていないのでしょう。


偏微分方程式を解くというより、解かりきった答えを、sf.solvePDE(..) で再現したにすぎない、数値解析結果です。でも、ここに到達するまでに 20 時間以上、この方程式と格闘しています。その過程で、この Hamilton-Jacobi 方程式が初期条件の設定しだいで直ぐに複素数値解に入り込むことを体験させられました。非線形な偏微分方程式では許される初期条件が大きく限られることを身をもって理解できました。解かりきった答えの再現になった理由を理解願います。

FDTD への適用: 二重スリットによる電磁波の干渉

電磁場数値解析の手法である FDTD(Finite Difference Time Domain) を itrSlvPDE(..) で解いてみます。FDTD は漸化式による偏微分方程式の数値解法とみなせます。solvePDE(..), itrSlvPDE(..) で扱えます。二重スリットの例の場合は itrSlvPDE(..) イタレータの方が変数スコープを同じにできるので少し簡潔に記述できます。

FDTD における漸化式

二重スリットによる干渉のための漸化式は、ここのJava ソースにあるものを使わせてもらいました。

FDTD の数値解

先の web page にある Java コードは干渉波形のアニメーション表示を行うコードですが、ここでは複数の波形データを HDD に落としておき, 後で pg.bat を使って 3D グラフとして表示させています。コードを単純化することを優先しています。

//@@
# coding=shift_jis
# 07.10.24 
# slit FDTD
import sf
import math

dt = 1.0
NNx = 201
NNy = 201

omega = math.pi/16.0
theta = 0.0
sqrt2v2 = 0.70710678

dctBoundary = {}    # 境界条件を設定する辞書
# とりあえず 30x30 行列の全部を 計算可能:True にする
for index in sf.mrng(NNx,NNy):
    dctBoundary[index] = True

# fncPDE(..) 内部でユーザー側で設定する。dctVctH の更新が (0,k) 位置について必要なため

for j in range(NNy):    # Ez 駆動ラインを 'i' にして特別扱いする
    dctBoundary[0,j] = 'i'

# スリット部分を指定します
for j in range(97,104):
    for k in range(0,64):
        dctBoundary[j,k] = False
    for k in range(80,120):
        dctBoundary[j,k] = False
    for k in range(136,NNy):
        dctBoundary[j,k] = False

dctVctH = {}            # 磁場ベクトル Hx,Hy を保持させる行列辞書
for index in sf.mrng(NNx+1,NNy+1):
    dctVctH[index] = [0,0]

countLoop = 0
lstDbg = [2,96] # to debug
import time
startTimeAt = time.time()
for dctRslt, index in sf.itrSlvPDE(
        dctBoundary, inIteration=251, blInterpolation = False ):
    i,j = index
    boundaryValue = dctBoundary[index]
    if index == tuple(lstDbg):  # to debug
        print "countLoop:", countLoop
    if (boundaryValue is True):
        #Ez[i][j] += sqrt2v2*((Hy[i+1][j]-Hy[i][j])-(Hx[i][j+1]-Hx[i][j]));
        dctRslt[index] += sqrt2v2*( (dctVctH[i+1,j][1]- dctVctH[i,j][1])
                                  - (dctVctH[i,j+1][0]- dctVctH[i,j][0]) )
    elif (boundaryValue is 'i'):
        if index == (0,0):
            theta += omega*dt
            if (theta<3.1416):
                # generate Ez at position [0,y]
                Ezt = math.sin(theta)
                for j in range(NNy):
                    dctRslt[0,j]=Ezt
    
            for i,j in sf.mrng(NNx,[1,NNy]):
                dctVctH[i,j][0] -= sqrt2v2*( dctRslt[i,j] - dctRslt[i,j-1])
            for i,j in sf.mrng([1,NNx],NNy):
                if (i,j) == tuple(lstDbg):  # to debug
                    pass
                dctVctH[i,j][1] += sqrt2v2*( dctRslt[i,j] - dctRslt[i-1,j])
    
            dctAt ={}
            for i,j in sf.mrng(NNx/5,NNy/5):
                dctAt[i,j] = dctRslt[5*i, 5*j]
    
            if countLoop % 50 == 0:
                sf.putSf(dctAt, "rslt"+str(countLoop) )
            countLoop += 1
print "Total executing time:",time.time() - startTimeAt
//@@@
//copy \#####.### temp.py /y

上のコードで自由端の線形外挿補間を禁止しています。blInterpolation = False のところです。これがないと数値解が発散します。FDTD が電場 E と磁場 E が交互に相手を変化させるもののためです。ここに自由端での線形補正を加えると数値解が発散します。

なお、上のコードでは、先の web page の Java コードにはあった端面での反射を防ぐコードを省略してあります。電磁波形がどのように伝わるのかを検討するときに無反射条件を追加するのは検討を複雑にするだけだからです。領域を広げれば無反射での挙動を見れるからです。そのため、上の NNx = 201,NNy = 201 で記述してあるコードでは、繰り返し回数を反射の影響のない 251 回までに限定しておく必要があります。

pg rslt100
index:(13, 0)  max: 1.01753
index:(9, 0)  min:-0.017216
pg rslt150
index:(20, 14)  max: 0.855047
index:(19, 0)  min:-0.770715
pg rslt200
index:(13, 20)  max: 0.356003
index:(12, 14)  min:-0.916944
pg rslt250
index:(6, 7)  max: 0.242679
index:(5, 14)  min:-0.916357


Wyle 方程式: ベクトル値偏微分方程式

ベクトル値の偏微分方程式であ Wyle 方程式を扱います。質量 0 のニュートリノが満たす方程式だと言われる下の方程式です。
  ∂ψ/∂t +   σx∂ψ/∂x +   σy ∂ψ/∂y +   σz ∂ψ/∂z == 0
ただし σx,σy,σz は Pauli 行列といわれる下の行列です
< 0,  1 >:σx
< 1,  0 >
sf "`σy"
< 0, -i >:σx
< i,  0 >

sf "`σz"
<  1,  0 >:σz
<  0, -1 >
ここで ψ(t,x) は二つの要素からなるベクトル値関数です。まずは一次元の Wyle 方程式を考えます。

∂ψ/∂t + σx∂ψ/∂x == 0 漸化式

∂ψ/∂t[j,k] ≒ (ψ[j,k+1] - ψ[j,k])/Δt ------------ (1 式)
∂ψ/∂x[j,k] ≒ (ψ[j+1,k] - ψ[j,k])/Δx ------------ (2a 式)
∂ψ/∂x[j,k] ≒ (ψ[j,k] - ψ[j-1,k])/Δx ------------ (2b 式)

0 == (1 式) + σx(2a 式)
  == (ψ[j,k+1] - ψ[j,k])/Δt + σx(ψ[j+1,k] - ψ[j,k])/Δx
∴
  ψ[j,k+1] == ψ[j,k] + σx(ψ[j+1,k] - ψ[j,k])Δt/Δx
  ψ[j,k] == ψ[j,k-1] + σx(ψ[j+1,k-1] - ψ[j,k-1])Δt/Δx --- (3a 式)

0 == (1 式) + (2b 式)
  == (ψ[j,k+1] - ψ[j,k])/Δt + σx(ψ[j,k] - ψ[j-1,k])/Δx

∴
  ψ[j,k] == ψ[j,k-1] + σx(ψ[j,k-1] - ψ[j-1,k-1])Δt/Δx --- (3b 式)

[(3a 式)+(3b 式)]/2 より
  ψ[j,k] == ψ[j,k-1] + 0.5σx(ψ[j+1,k-1] - ψ[j-1,k-1])Δt/Δx --- (3 式)

∂ψ/∂t + σx∂ψ/∂x == 0 の数値解 1

//@@
# 07.11.07
import sf
N = 200 # t
M = 60 # x
sigmaX = sf.sc.array([ [0,1]
                      ,[1,0] ])
deltaT = 0.1
#deltaT = 1.0
deltaX = 1.0
dctBoundary = {}
for index in sf.mrng(N,M):  dctBoundary[index] = True

pi = sf.math.pi
for i in range(M):
    r = sf.math.exp(-((1.0*i-M/2)/5 )**2)
    dctBoundary[0,i] = sf.sc.array([r,0])


lstDbg = [2,1]
inCount = 0
beforeTimeIndex = 0
for dctRslt, index in sf.itrSlvPDE(dctBoundary, dfltValue = sf.sc.array([0.0,0.0])
#         , inIteration =  1 ):
        , inIteration = 10 ):
    if tuple(lstDbg) == index:
        print "loopCount:",inCount
        inCount+=1
    j,k = index
    dctRslt[j,k] = dctRslt[j-1,k] + 0.5*sf.sc.dot(sigmaX,
        (dctRslt[j-1,k+1] - dctRslt[j-1,k-1])*deltaT/deltaX )

# display result
clAt = sf.ClCplxColor()
dctAt ={}
#K=10
K=1
L=4
#L=10
for index in sf.mrng(N,M):
    j,k=index
    if (k % K != 0) or (j%L != 0):
        continue
    
    z = dctRslt[index][0]+ dctRslt[index][1]*1j
    if z == 0+0j:
        phase = 0+0j
    else:
        phase = 0.999*z/abs(z)

    dctAt[k//K,j//L,] = sf.vs.mag( dctRslt[index] )

sf.make3dRowCol(dctAt)
sf.putSf(dctRslt)
//@@@
//copy \#####.### temp2.py /y

上のプログラムで deltaT = 0.1, deltaX = 1.0 と時間と位置で 10 倍の比としています。(これを 1:1 の比にすると計算結果は発散します。また波が端に到達して反射するようになっても計算結果が発散するのでご注意ください。)その結果として時間と位置の分割非を 200:60 としたので、表示するときは時間軸側を四つごとに抜き出して表示させています。下のような結果となります。

時刻 0 でガウス分布していたニュートリノが左右に分裂して光速度で飛んで行っています。左右に分かれた後の波束は、拡散することなく同じ形を保っています。

t,x 平面上でのベクトル関数 ψ(t,x) の時間変化を一つの 3D グラフにするため sf.vs.mag( dctRslt[index] ) と ψ(t,x):dctRslt[index] ベクトル値の絶対値:norm を表示させるようにしています。sf.putSf(dctRslt) によって計算結果が _rs.val という sf 変数ファイルに落としてあるので、時刻 t = 0.1*100 の位置での ψ(0.1*100,x) ベクトル値がどう分布しているのかを、下のように sf を使って数値を直接見ることができます。

N=60,<>
< -1.23001e-006, -1.23001e-006 >
<  3.51873e-006,  3.51873e-006 >
<  8.20815e-006,  8.20815e-006 >
<  2.46361e-005,  2.46361e-005 >
<  6.80837e-005,  6.80837e-005 >
<   0.000180167,   0.000180167 >
<   0.000454054,   0.000454054 >
<    0.00108905,    0.00108905 >
<     0.0024828,     0.0024828 >
<    0.00537315,    0.00537315 >
<     0.0110231,     0.0110231 >
<     0.0214053,     0.0214053 >
<     0.0392807,     0.0392807 >
<     0.0680021,     0.0680021 >
<      0.110847,      0.110847 >
<      0.169779,      0.169779 >
<      0.243783,      0.243783 >
<      0.327317,      0.327317 >
<      0.409757,      0.409757 >
<      0.476701,      0.476701 >
<      0.513416,      0.513416 >
<      0.509609,      0.509609 >
<      0.463636,      0.463637 >
<      0.384003,      0.384002 >
<      0.287002,       0.28699 >
<      0.191256,      0.191222 >
<      0.111631,      0.111629 >
<     0.0553369,     0.0555941 >
<     0.0217501,     0.0226481 >
<    0.00548843,    0.00688825 >
<   0.000854666,             0 >
<    0.00548843,   -0.00688825 >
<     0.0217501,    -0.0226481 >
<     0.0553369,    -0.0555941 >
<      0.111631,     -0.111629 >
<      0.191256,     -0.191222 >
<      0.287002,      -0.28699 >
<      0.384003,     -0.384002 >
<      0.463636,     -0.463637 >
<      0.509609,     -0.509609 >
<      0.513416,     -0.513416 >
<      0.476701,     -0.476701 >
<      0.409757,     -0.409757 >
<      0.327317,     -0.327317 >
<      0.243783,     -0.243783 >
<      0.169779,     -0.169779 >
<      0.110847,     -0.110847 >
<     0.0680021,    -0.0680021 >
<     0.0392807,    -0.0392807 >
<     0.0214053,    -0.0214053 >
<     0.0110231,    -0.0110231 >
<    0.00537315,   -0.00537315 >
<     0.0024828,    -0.0024828 >
<    0.00108904,   -0.00108904 >
<   0.000454079,  -0.000454079 >
<   0.000180082,  -0.000180082 >
<  6.83517e-005, -6.83517e-005 >
<  2.38295e-005, -2.38295e-005 >
<  1.05205e-005, -1.05205e-005 >
< -2.78856e-006,  2.78856e-006 >

二つの山ができていることが、上の数値列から分ります。またベクトル値二成分の比率も概ね一定のようです。

∂ψ/∂t + σx∂ψ/∂x == 0 の数値解 2

上の数値解のグラフ表示では、ベクトル値関数 ψ(x,t)の二成分ベクトルがどのように変化しているのか直接読み取れません。sf で時刻を指定して数値列を抜き出していたのでは、全体でのベクトル二成分の変化の様子は読み取れません。_rs.val 全部の数字を打ち出しても、数字の塊が出てくるだけであり、全体のベクトル成分の様子を掴めません。

幸いにも上の ψ(x,t) の二成分ベクトルは実数値です。二つの実数値は複素数とみなすことで一つの数値にできます。さらに幸いなことに、複素数値の平面分布での位相変化を kkRGB と名づけた RGB 三色の変化で表示させる機能を sf.py に用意してあります。それを 3D 表示させる makeFacesWithRGB(..) 関数も sf.py に用意してあります。それを使って一次元 Wyle 方程式の数値解のベクトル成分変化を表示させて見ましょう。今度の初期状態は位相変化をもつ矩形とします。

//@@
# 07.11.07
import sf
N = 200 # t
M = 60 # x
sigmaX = sf.sc.array([ [0,1]
                      ,[1,0] ])
deltaT = 0.1
#deltaT = 1.0
deltaX = 1.0
dctBoundary = {}
for index in sf.mrng(N,M):  dctBoundary[index] = True

pi = sf.math.pi
"""'
for i in range(M):
    r = sf.math.exp(-((1.0*i-M/2)/5 )**2)
    dctBoundary[0,i] = sf.sc.array([r,0])
'"""
for i in range(10):
    dctBoundary[0,25+i] = sf.sc.array([sf.math.cos(2*pi*i/10), sf.math.sin(2*pi*i/10)])
    #dctBoundary[0,i] = sf.sc.array([1,0])


lstDbg = [2,1]
inCount = 0
beforeTimeIndex = 0
for dctRslt, index in sf.itrSlvPDE(dctBoundary, dfltValue = sf.sc.array([0.0,0.0])
#         , inIteration =  1 ):
        , inIteration = 10 ):
    if tuple(lstDbg) == index:
        print "loopCount:",inCount
        inCount+=1
    j,k = index
    dctRslt[j,k] = dctRslt[j-1,k] + 0.5*sf.sc.dot(sigmaX,
        (dctRslt[j-1,k+1] - dctRslt[j-1,k-1])*deltaT/deltaX )

# display result
clAt = sf.ClCplxColor()
dctAt ={}
#K=10
K=1
L=4
#L=10
for index in sf.mrng(N,M):
    j,k=index
    if (k % K != 0) or (j%L != 0):
        continue
    
    z = dctRslt[index][0]+ dctRslt[index][1]*1j
    if z == 0+0j:
        phase = 0+0j
    else:
        phase = 0.999*z/abs(z)

    dctAt[k//K,j//L,] = ( sf.sc.array((k//K,j//L,10*abs(z)))-[M/(2*K),N/(2*L),0], clAt.GetFltColor(phase) )
    #dctAt[k//K,j//L] = ( sf.sc.array((k//K,j//L,10*abs(z)))-[N/(2*K),M/(2*L),0], clAt.GetFltColor(phase) )

sf.makeFacesWithRGB(dctAt, blMesh = True, blDirection = False)
sf.putSf(dctRslt)
//@@@
//copy \#####.### temp2.py /y

計算結果の kkRGB は下のようになります。

特定の位相:二成分ベクトルの特定の比成分が左右に分かれて飛んでいく成分となります。左右に分かれた後の二成分ベクトルの比は一定であることが判ります。