［境界要素法プログラムを設計する（計算部と出力部の実装）］ [ネコ騙し数学]
［境界要素法プログラムを設計する（計算部と出力部の実装）］
１．実行者は要素だ。そして常にトップダウンだ！
入力部の実装によって、計算に必要な全ての情報はセットされました（されたはず(^^;)）。次にやる事は、それらを利用して具体的な計算ルーティンを書く事です。計算部（ビジネスルール層）は、アプリケーションの中で最も動的な部分で、書いてて最も面白い部分でもありますが、バグの温床にもなります。最も動的だからです。
人間は動的対象を見るとすぐに幻惑され、目算を見失います。なので動的過程を動かす、メタな静的構造を保持しておくのが安全です。それは多くの場合、動的過程の制御部分ですが、これとトップダウン思考を重ねると、制御部分は「実行者」になります。「実行者」は「要素」でした。
[計算部]の冒頭にこうあります。
REM 境界要素番号kで境界積分を行うLoop
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
要するにLoopを作りますが、設計の基本は常トップダウンです。よってLoopは、最外殻から書き始めるのが基本です。そうして動的過程を動かすメタな静的構造を確定させると、そこに「実行者」がすんなり納まります。要するに、境界要素を走るLoopを最初に書くべきだ、という事です。
REM 境界要素番号kで境界積分を行うLoop
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
Next k
ところで上記を見て、何か気づきませんか？。目標とする方程式系は、
そこで慌てて[メモ欄]を見返すと、定義してない！。[メモ欄]と[宣言部]は一対一に対応してるので、要はAを宣言してないって事です。そこでさらに慌てて、[メモ欄]のz(i)に関する部分の後に、
REM A(i, j)：全体係数マトリックス、i,j=1～2*BN_Count
と書き加え同時に、[宣言部]の同じ場所に、
Dim A(1000, 1000) REM 全体係数マトリックス、i,j=1～2*BN_Count
を追加します。i,j=1～2*BN_Countなのは、z(i)，i=1～2*BN_Countだからです。
実務でもこんな事は良くあるのだ！。じっさい忘れてたし(^^;)。でもトップダウン方式だから被害は最小限と思えませんか？。これがボトムアップだったら、既に変数Aは別の目的で使われてて2重宣言が起こり、いや2度ある事は3度ある・・・3重宣言が起こり、3つあるAの内の2つをA1とA2に書き変える事にして、沢山の行のAを一つ一つ判断しながらA1やA2に書き変え、そこでまた判断ミスしてそれがまたバグの温床となり、・・・で訳わかんなくなってバグの連鎖反応が、・・・などと屁理屈をこねるのであった(^^;)。
全体係数マトリックスAを導入したおかげで、さらにもう一つの大事に気づきます。思い起こせば境界積分の値は、特異点の位置によって変わるのでした。マトリックスAは境界方程式に対応するので、いま考える特異点は境界節点のどれかと一致します。その番号（特異点番号）をiとした時、iはAに含まれるBとHの行番号なのでした。つまり、
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
Next k
Next i
としなければならない訳です。まぁ～、こういう事も良くありますって・・・(^^;)。
次にもう少し「境界要素番号kで境界積分を行うLoop」を具体化します。入力部の情報をいかに使うべきかを、自分にわからせるためにです。
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
x0 = BN(i, 1) REM 節点iの(x，y)座標
y0 = BN(i, 2)
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
j1 = EN_B(k，1) REM 要素kの始端の節点番号
j2 = EN_B(k，2) REM 要素kの終端の節点番号
x1 = BN(j1, 1) REM 節点j1の(x，y)座標
y1 = BN(j1, 2)
x2 = BN(j2, 1) REM 節点j2の(x，y)座標
y2 = BN(j2, 2)
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
jj1 = BN_Z(j1，1) REM ψj1の解ﾍﾞｸﾄﾙzの中での位置
jjj1 = BN_Z(j1，2) REM qj1の解ﾍﾞｸﾄﾙzの中での位置
jj2 = BN_Z(j2，1) REM ψj2の解ﾍﾞｸﾄﾙzの中での位置
jjj2 = BN_Z(j2，2) REM qj2の解ﾍﾞｸﾄﾙzの中での位置
A(i, jj1) = A(i, jj1) - B1
A(i, jjj1) = A(i, jjj1) + H1
A(i, jj2) = A(i, jj2) - B2
A(i, jjj2) = A(i, jjj2) + H2
Next k
Next i
赤字の部分はいま追加したコメントで、後で考える必要のある部分です。
「A(i, jj1) = A(i, jj1) - B1」などの部分は、次の理由によります。一つの特異点iについて境界積分は具体的には、次の形になりました。
(ξi，ηi)を番号（節点番号）iの特異点だとして、
A(i, jj1) = A(i, jj1) - B1
A(i, jj2) = A(i, jj2) - B2
などとしとけば良いじゃん！、という話です。こうしとけば、「For k = 1 to B_Count」のLoopを回しただけで、自然に要素kとk+1に関する「重ね合わせ（足し算）」が起きるからです。そのためにこそ、要素－節点番号対応表EN_B(k，2)があるんですよ。
※ 宣言時にA(i, j)は、自動で全ての値が0に初期化されるのが普通です。
有限要素法などだと少なくとも2次元なので、要素－節点番号対応表がないとプログラムすら書けません。要素－節点番号対応表は、2次元領域を有限要素分割して初めて決まります。そういう訳でメッシュ切り（有限要素分割）と要素－節点番号対応表の生成は、FEMの中で最も手のかかる作業なんです、現実には(^^;)。
とりあえず境界方程式の核心部分は、これで生成できるはずです。次へ進みましょう。次には、
REM 領域要素番号kで領域積分を行うLoop
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
と書いてあります。念のため次の事を注意します。
これから最後のREM文に書いてある既知ﾍﾞｸﾄﾙwをセットする事になりますが、wの値は解ﾍﾞｸﾄﾙ（未知ﾍﾞｸﾄﾙ）zへセットします。これは連立方程式プログラムでは普通に行われる事で、係数マトリックスAと既知ﾍﾞｸﾄﾙとしてのzを連立方程式の解法部（たいていサブルーティン）へ送ると、zに解がセットされて返って来る仕掛けになるからです。こうすると既知ﾍﾞｸﾄﾙ分のメモリを節約できます。
ここも基本は同じです。
REM 領域要素番号kで領域積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
h = 0
For k = 1 to R_Count
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 特異点iと要素節点情報から積分値w0を計算
h = h + w0
Next k
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
z(i) = h
Next i
領域積分項wは、ラプラス方程式Δψ＝gの既知関数g(x，y)に対応するものですが、今回g(x，y)は具体的に与えてないので、領域積分についてはこれ以上具体化しません。プログラムテストのところで、いくつか具体的例を出す予定です。
これで離散化された境界方程式はセットされました（されたはず(^^;)）。残りは境界条件のセットです。
REM 境界節点番号jで境界条件位置を指定するLoop
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
・・・なんですが、既に最外殻Loopは、全体係数マトリックスAの行番号iに関するものなのがわかっています。境界条件は境界方程式に対する補助条件なので、Aの中で境界方程式の行並びの後に、単純に並べれば良いだけです。境界方程式の行並びは特異点番号に一致し、それは特異点が境界節点のどれかと一致する事を表す並びでもあったので、境界方程式の行並びの行番号はi = 1～BN_Countでした。境界条件の行並びはi = BN_Count + 1から始まります。境界条件は1個の境界節点に必ず1個あります。従って、
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
Next i
いまLoop内のREM文の中の節点番号jは、さっきの説明から、j＝i－BN_Countです。そこに注意すると、
※「十進BASIC 最新版 概要」によると、IFブロックが使えます。
http://hp.vector.co.jp/authors/VA008683/basic02.htm
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 行番号から節点番号を特定する
j = i - BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
If BN_BP(j，1) = 1 Then
REM ψjが既知が場合
jj = BN_Z(j，1) REM ψjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM ψjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，1) REM ψjの値
ElseIf BN_BP(j，2) = 1 Then
REM qjが既知が場合
jj = BN_Z(j，2) REM qjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM qjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，2) REM qjの値
End If
Next i
IFブロックを書く時もトップダウンの作法に従い、まずIf～ElseIf～End Ifのブロックの形を決めてから、必要なら説明のためのREM分をブロック内に書き込み、最後に具体的なコードを書く事をお奨めします。
最後に、全体係数マトリックスAと既知ﾍﾞｸﾄﾙとしてのzを連立方程式の解法部へ送ります。
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ﾍﾞｸﾄﾙzを求める
これも後で考えるべきものとして、赤字化しておきます。
[出力部]は短いので、それもここでやってしまいましょう。
REM BN_Z(j，2) に従いz(j)から値を読み出し、出力
For j = 1 to BN_Count
jj = BN_Z(j，1) REM zの中でのψjの位置
REM z(jj)をψjの出力ﾙｰﾃｨﾝへ渡す
Next i
For j = 1 to BN_Count
jj = BN_Z(j，2) REM zの中でのqjの位置
REM z(jj)をqjの出力ﾙｰﾃｨﾝへ渡す
Next i
赤字の部分は、サブルーティン化する予定です。でもその前に、本当にサブルーティン化する必要があるのか？を検討して、損はありません。
REM [メモ欄]
REM B_Count：境界要素数
REM EN_B(k，2)：境界要素－節点対応表，k=1～B_Count
REM (k，1)：要素kの始端の節点No.，(k，2)：終端の節点No.
REM R_Count：領域要素数
REM EN_R(k，4)：領域要素－節点対応表，k=1～R_Count
REM (k，j)，j=1～4：jに従って要素kの要素節点No.を左回りに格納
REM BN_Count：境界節点数
REM BN(j，2)：境界節点座標，j=1～BN_Count，(j，1)：節点jのx座標，(j，2)：y座標
REM RN_Count：領域節点数
REM RN(j，2)：領域節点座標，j=1～RN_Count，(j，1)：節点jのx座標，(j，2)：y座標
REM BN_BP(j，2)：節点の境界条件の有無，j=1～BN_Count
REM (j，1)：節点jのψjが既知なら1、そうでないなら0
REM (j，2)：節点jのqjが既知なら1、そうでないなら0
REM BN_BV(j，2)：境界条件の値，j=1～BN_Count
REM (j，1)：節点jのψjに与える値
REM (j，2)：節点jのqjに与える値
REM z(i):未知ベクトル，i=1～2*BN_Count
REM BN_Z(j，2)：節点－未知数対応表，j=1～BN_Count
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)
REM A(i, j)：全体係数マトリックス、i,j=1～2*BN_Count
REM [宣言部]
REM 配列寸法の上限は1000
REM B_Count：境界要素数
REM R_Count：領域要素数
Dim EN_B(1000，2) REM 境界要素－節点対応表，k=1～B_Count
Dim EN_R(1000，4) REM 領域要素－節点対応表，k=1～R_Count
REM BN_Count：境界節点数
REM RN_Count：領域節点数
Dim BN(1000，2) REM 境界節点座標，j=1～BN_Count
Dim RN(1000，2) REM 領域節点座標，j=1～RN_Count
Dim BN_BP(1000，2) REM 節点の境界条件の有無，j=1～BN_Count
Dim BN_BV(1000，2) REM 境界条件の値，j=1～BN_Count
Dim BN_Z(1000，2) REM 節点－未知数対応表，j=1～BN_Count
Dim z(1000) REM 未知ベクトル，i=1～2*BN_Count
Dim A(1000, 1000) REM 全体係数マトリックス、i,j=1～2*BN_Count
REM [入力部]
REM境界要素数B_Countと領域要素数R_Countは、事前に与えられると仮定する。
REM 境界要素と領域要素データは、要素ごとに節点番号が与えられると仮定する。
REM 与えられる節点番号は、左回りに順序付けられたものが与えられると仮定する。
REM 境界要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、境界要素－節点対応表EN_B(k，2)をセット
For k = 1 to B_Count
EN_B(k，1) = ？
EN_B(k，2) = ？
Next k
REM 領域要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、領域要素－節点対応表EN_R(k，4)をセット
For k = 1 to R_Count
For j = 1 to 4
EN_R(k，j) = ？
Next j
Next k
REM境界節点数BN_Countと領域節点数RN_Countは、事前に与えられると仮定する。
REM 境界節点を節点番号jで読み込むLoop
REM Loop内で境界節点座標を読み、BN(j，2)をセット
For j = 1 to BN_Count
BN(j，1) = ？
BN(j，2) = ？
Next j
REM 領域節点を節点番号jで読み込むLoop
REM Loop内で領域節点座標を読み、RN(j，2)をセット
For j = 1 to RN_Count
RN(j，1) = ？
RN(j，2) = ？
Next j
REM 境界条件の場所を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の有無を読み、BN_BP(j，2)をセット
For j = 1 to BN_Count
BN_BP(j，1) = ？
BN_BP(j，2) = ？
Next j
REM 境界条件の値を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の値を読み、BN_BV(j，2)をセット
For j = 1 to BN_Count
BN_BV(j，1) = ？
BN_BV(j，2) = ？
Next j
REM 節点番号jのLoopで、節点－未知数対応表BN_Z(j，2)をセット
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 対応は、節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)
For j = 1 to BN_Count
BN_Z(j，1) = j
BN_Z(j，2) = BN_Count + j
Next j
REM [計算部]
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
x0 = BN(i, 1) REM 節点iの(x，y)座標
y0 = BN(i, 2)
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
j1 = EN_B(k，1) REM 要素kの始端の節点番号
j2 = EN_B(k，2) REM 要素kの終端の節点番号
x1 = BN(j1, 1) REM 節点j1の(x，y)座標
y1 = BN(j1, 2)
x2 = BN(j2, 1) REM 節点j2の(x，y)座標
y2 = BN(j2, 2)
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
jj1 = BN_Z(j1，1) REM ψj1の解ﾍﾞｸﾄﾙzの中での位置
jjj1 = BN_Z(j1，2) REM qj1の解ﾍﾞｸﾄﾙzの中での位置
jj2 = BN_Z(j2，1) REM ψj2の解ﾍﾞｸﾄﾙzの中での位置
jjj2 = BN_Z(j2，2) REM qj2の解ﾍﾞｸﾄﾙzの中での位置
A(i, jj1) = A(i, jj1) - B1
A(i, jjj1) = A(i, jjj1) + H1
A(i, jj2) = A(i, jj2) - B2
A(i, jjj2) = A(i, jjj2) + H2
Next k
Next i
REM 領域要素番号kで領域積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
h = 0
For k = 1 to R_Count
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 特異点iと要素節点情報から積分値w0を計算
h = h + w0
Next k
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
z(i) = h
Next i
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 行番号から節点番号を特定する
j = i - BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
If BN_BP(j，1) = 1 Then
REM ψjが既知が場合
jj = BN_Z(j，1) REM ψjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM ψjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，1) REM ψjの値
ElseIf BN_BP(j，2) = 1 Then
REM qjが既知が場合
jj = BN_Z(j，2) REM qjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM qjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，2) REM qjの値
End If
Next i
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ﾍﾞｸﾄﾙzを求める
REM [出力部]
REM BN_Z(j，2) に従いz(j)から値を読み出し、出力
For j = 1 to BN_Count
jj = BN_Z(j，1) REM zの中でのψjの位置
REM z(jj)をψjの出力ﾙｰﾃｨﾝへ渡す
Next i
For j = 1 to BN_Count
jj = BN_Z(j，2) REM zの中でのqjの位置
REM z(jj)をqjの出力ﾙｰﾃｨﾝへ渡す
Next i
１．実行者は要素だ。そして常にトップダウンだ！
入力部の実装によって、計算に必要な全ての情報はセットされました（されたはず(^^;)）。次にやる事は、それらを利用して具体的な計算ルーティンを書く事です。計算部（ビジネスルール層）は、アプリケーションの中で最も動的な部分で、書いてて最も面白い部分でもありますが、バグの温床にもなります。最も動的だからです。
人間は動的対象を見るとすぐに幻惑され、目算を見失います。なので動的過程を動かす、メタな静的構造を保持しておくのが安全です。それは多くの場合、動的過程の制御部分ですが、これとトップダウン思考を重ねると、制御部分は「実行者」になります。「実行者」は「要素」でした。
[計算部]の冒頭にこうあります。
REM 境界要素番号kで境界積分を行うLoop
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
要するにLoopを作りますが、設計の基本は常トップダウンです。よってLoopは、最外殻から書き始めるのが基本です。そうして動的過程を動かすメタな静的構造を確定させると、そこに「実行者」がすんなり納まります。要するに、境界要素を走るLoopを最初に書くべきだ、という事です。
REM 境界要素番号kで境界積分を行うLoop
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
Next k
ところで上記を見て、何か気づきませんか？。目標とする方程式系は、
そこで慌てて[メモ欄]を見返すと、定義してない！。[メモ欄]と[宣言部]は一対一に対応してるので、要はAを宣言してないって事です。そこでさらに慌てて、[メモ欄]のz(i)に関する部分の後に、
REM A(i, j)：全体係数マトリックス、i,j=1～2*BN_Count
と書き加え同時に、[宣言部]の同じ場所に、
Dim A(1000, 1000) REM 全体係数マトリックス、i,j=1～2*BN_Count
を追加します。i,j=1～2*BN_Countなのは、z(i)，i=1～2*BN_Countだからです。
実務でもこんな事は良くあるのだ！。じっさい忘れてたし(^^;)。でもトップダウン方式だから被害は最小限と思えませんか？。これがボトムアップだったら、既に変数Aは別の目的で使われてて2重宣言が起こり、いや2度ある事は3度ある・・・3重宣言が起こり、3つあるAの内の2つをA1とA2に書き変える事にして、沢山の行のAを一つ一つ判断しながらA1やA2に書き変え、そこでまた判断ミスしてそれがまたバグの温床となり、・・・で訳わかんなくなってバグの連鎖反応が、・・・などと屁理屈をこねるのであった(^^;)。
全体係数マトリックスAを導入したおかげで、さらにもう一つの大事に気づきます。思い起こせば境界積分の値は、特異点の位置によって変わるのでした。マトリックスAは境界方程式に対応するので、いま考える特異点は境界節点のどれかと一致します。その番号（特異点番号）をiとした時、iはAに含まれるBとHの行番号なのでした。つまり、
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
Next k
Next i
としなければならない訳です。まぁ～、こういう事も良くありますって・・・(^^;)。
次にもう少し「境界要素番号kで境界積分を行うLoop」を具体化します。入力部の情報をいかに使うべきかを、自分にわからせるためにです。
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
x0 = BN(i, 1) REM 節点iの(x，y)座標
y0 = BN(i, 2)
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
j1 = EN_B(k，1) REM 要素kの始端の節点番号
j2 = EN_B(k，2) REM 要素kの終端の節点番号
x1 = BN(j1, 1) REM 節点j1の(x，y)座標
y1 = BN(j1, 2)
x2 = BN(j2, 1) REM 節点j2の(x，y)座標
y2 = BN(j2, 2)
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
jj1 = BN_Z(j1，1) REM ψj1の解ﾍﾞｸﾄﾙzの中での位置
jjj1 = BN_Z(j1，2) REM qj1の解ﾍﾞｸﾄﾙzの中での位置
jj2 = BN_Z(j2，1) REM ψj2の解ﾍﾞｸﾄﾙzの中での位置
jjj2 = BN_Z(j2，2) REM qj2の解ﾍﾞｸﾄﾙzの中での位置
A(i, jj1) = A(i, jj1) - B1
A(i, jjj1) = A(i, jjj1) + H1
A(i, jj2) = A(i, jj2) - B2
A(i, jjj2) = A(i, jjj2) + H2
Next k
Next i
赤字の部分はいま追加したコメントで、後で考える必要のある部分です。
「A(i, jj1) = A(i, jj1) - B1」などの部分は、次の理由によります。一つの特異点iについて境界積分は具体的には、次の形になりました。
(ξi，ηi)を番号（節点番号）iの特異点だとして、
A(i, jj1) = A(i, jj1) - B1
A(i, jj2) = A(i, jj2) - B2
などとしとけば良いじゃん！、という話です。こうしとけば、「For k = 1 to B_Count」のLoopを回しただけで、自然に要素kとk+1に関する「重ね合わせ（足し算）」が起きるからです。そのためにこそ、要素－節点番号対応表EN_B(k，2)があるんですよ。
※ 宣言時にA(i, j)は、自動で全ての値が0に初期化されるのが普通です。
有限要素法などだと少なくとも2次元なので、要素－節点番号対応表がないとプログラムすら書けません。要素－節点番号対応表は、2次元領域を有限要素分割して初めて決まります。そういう訳でメッシュ切り（有限要素分割）と要素－節点番号対応表の生成は、FEMの中で最も手のかかる作業なんです、現実には(^^;)。
とりあえず境界方程式の核心部分は、これで生成できるはずです。次へ進みましょう。次には、
REM 領域要素番号kで領域積分を行うLoop
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
と書いてあります。念のため次の事を注意します。
これから最後のREM文に書いてある既知ﾍﾞｸﾄﾙwをセットする事になりますが、wの値は解ﾍﾞｸﾄﾙ（未知ﾍﾞｸﾄﾙ）zへセットします。これは連立方程式プログラムでは普通に行われる事で、係数マトリックスAと既知ﾍﾞｸﾄﾙとしてのzを連立方程式の解法部（たいていサブルーティン）へ送ると、zに解がセットされて返って来る仕掛けになるからです。こうすると既知ﾍﾞｸﾄﾙ分のメモリを節約できます。
ここも基本は同じです。
REM 領域要素番号kで領域積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
h = 0
For k = 1 to R_Count
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 特異点iと要素節点情報から積分値w0を計算
h = h + w0
Next k
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
z(i) = h
Next i
領域積分項wは、ラプラス方程式Δψ＝gの既知関数g(x，y)に対応するものですが、今回g(x，y)は具体的に与えてないので、領域積分についてはこれ以上具体化しません。プログラムテストのところで、いくつか具体的例を出す予定です。
これで離散化された境界方程式はセットされました（されたはず(^^;)）。残りは境界条件のセットです。
REM 境界節点番号jで境界条件位置を指定するLoop
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
・・・なんですが、既に最外殻Loopは、全体係数マトリックスAの行番号iに関するものなのがわかっています。境界条件は境界方程式に対する補助条件なので、Aの中で境界方程式の行並びの後に、単純に並べれば良いだけです。境界方程式の行並びは特異点番号に一致し、それは特異点が境界節点のどれかと一致する事を表す並びでもあったので、境界方程式の行並びの行番号はi = 1～BN_Countでした。境界条件の行並びはi = BN_Count + 1から始まります。境界条件は1個の境界節点に必ず1個あります。従って、
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
Next i
いまLoop内のREM文の中の節点番号jは、さっきの説明から、j＝i－BN_Countです。そこに注意すると、
※「十進BASIC 最新版 概要」によると、IFブロックが使えます。
http://hp.vector.co.jp/authors/VA008683/basic02.htm
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 行番号から節点番号を特定する
j = i - BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
If BN_BP(j，1) = 1 Then
REM ψjが既知が場合
jj = BN_Z(j，1) REM ψjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM ψjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，1) REM ψjの値
ElseIf BN_BP(j，2) = 1 Then
REM qjが既知が場合
jj = BN_Z(j，2) REM qjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM qjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，2) REM qjの値
End If
Next i
IFブロックを書く時もトップダウンの作法に従い、まずIf～ElseIf～End Ifのブロックの形を決めてから、必要なら説明のためのREM分をブロック内に書き込み、最後に具体的なコードを書く事をお奨めします。
最後に、全体係数マトリックスAと既知ﾍﾞｸﾄﾙとしてのzを連立方程式の解法部へ送ります。
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ﾍﾞｸﾄﾙzを求める
これも後で考えるべきものとして、赤字化しておきます。
[出力部]は短いので、それもここでやってしまいましょう。
REM BN_Z(j，2) に従いz(j)から値を読み出し、出力
For j = 1 to BN_Count
jj = BN_Z(j，1) REM zの中でのψjの位置
REM z(jj)をψjの出力ﾙｰﾃｨﾝへ渡す
Next i
For j = 1 to BN_Count
jj = BN_Z(j，2) REM zの中でのqjの位置
REM z(jj)をqjの出力ﾙｰﾃｨﾝへ渡す
Next i
赤字の部分は、サブルーティン化する予定です。でもその前に、本当にサブルーティン化する必要があるのか？を検討して、損はありません。
REM [メモ欄]
REM B_Count：境界要素数
REM EN_B(k，2)：境界要素－節点対応表，k=1～B_Count
REM (k，1)：要素kの始端の節点No.，(k，2)：終端の節点No.
REM R_Count：領域要素数
REM EN_R(k，4)：領域要素－節点対応表，k=1～R_Count
REM (k，j)，j=1～4：jに従って要素kの要素節点No.を左回りに格納
REM BN_Count：境界節点数
REM BN(j，2)：境界節点座標，j=1～BN_Count，(j，1)：節点jのx座標，(j，2)：y座標
REM RN_Count：領域節点数
REM RN(j，2)：領域節点座標，j=1～RN_Count，(j，1)：節点jのx座標，(j，2)：y座標
REM BN_BP(j，2)：節点の境界条件の有無，j=1～BN_Count
REM (j，1)：節点jのψjが既知なら1、そうでないなら0
REM (j，2)：節点jのqjが既知なら1、そうでないなら0
REM BN_BV(j，2)：境界条件の値，j=1～BN_Count
REM (j，1)：節点jのψjに与える値
REM (j，2)：節点jのqjに与える値
REM z(i):未知ベクトル，i=1～2*BN_Count
REM BN_Z(j，2)：節点－未知数対応表，j=1～BN_Count
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)
REM A(i, j)：全体係数マトリックス、i,j=1～2*BN_Count
REM [宣言部]
REM 配列寸法の上限は1000
REM B_Count：境界要素数
REM R_Count：領域要素数
Dim EN_B(1000，2) REM 境界要素－節点対応表，k=1～B_Count
Dim EN_R(1000，4) REM 領域要素－節点対応表，k=1～R_Count
REM BN_Count：境界節点数
REM RN_Count：領域節点数
Dim BN(1000，2) REM 境界節点座標，j=1～BN_Count
Dim RN(1000，2) REM 領域節点座標，j=1～RN_Count
Dim BN_BP(1000，2) REM 節点の境界条件の有無，j=1～BN_Count
Dim BN_BV(1000，2) REM 境界条件の値，j=1～BN_Count
Dim BN_Z(1000，2) REM 節点－未知数対応表，j=1～BN_Count
Dim z(1000) REM 未知ベクトル，i=1～2*BN_Count
Dim A(1000, 1000) REM 全体係数マトリックス、i,j=1～2*BN_Count
REM [入力部]
REM境界要素数B_Countと領域要素数R_Countは、事前に与えられると仮定する。
REM 境界要素と領域要素データは、要素ごとに節点番号が与えられると仮定する。
REM 与えられる節点番号は、左回りに順序付けられたものが与えられると仮定する。
REM 境界要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、境界要素－節点対応表EN_B(k，2)をセット
For k = 1 to B_Count
EN_B(k，1) = ？
EN_B(k，2) = ？
Next k
REM 領域要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、領域要素－節点対応表EN_R(k，4)をセット
For k = 1 to R_Count
For j = 1 to 4
EN_R(k，j) = ？
Next j
Next k
REM境界節点数BN_Countと領域節点数RN_Countは、事前に与えられると仮定する。
REM 境界節点を節点番号jで読み込むLoop
REM Loop内で境界節点座標を読み、BN(j，2)をセット
For j = 1 to BN_Count
BN(j，1) = ？
BN(j，2) = ？
Next j
REM 領域節点を節点番号jで読み込むLoop
REM Loop内で領域節点座標を読み、RN(j，2)をセット
For j = 1 to RN_Count
RN(j，1) = ？
RN(j，2) = ？
Next j
REM 境界条件の場所を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の有無を読み、BN_BP(j，2)をセット
For j = 1 to BN_Count
BN_BP(j，1) = ？
BN_BP(j，2) = ？
Next j
REM 境界条件の値を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の値を読み、BN_BV(j，2)をセット
For j = 1 to BN_Count
BN_BV(j，1) = ？
BN_BV(j，2) = ？
Next j
REM 節点番号jのLoopで、節点－未知数対応表BN_Z(j，2)をセット
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 対応は、節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)
For j = 1 to BN_Count
BN_Z(j，1) = j
BN_Z(j，2) = BN_Count + j
Next j
REM [計算部]
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
x0 = BN(i, 1) REM 節点iの(x，y)座標
y0 = BN(i, 2)
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
j1 = EN_B(k，1) REM 要素kの始端の節点番号
j2 = EN_B(k，2) REM 要素kの終端の節点番号
x1 = BN(j1, 1) REM 節点j1の(x，y)座標
y1 = BN(j1, 2)
x2 = BN(j2, 1) REM 節点j2の(x，y)座標
y2 = BN(j2, 2)
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
jj1 = BN_Z(j1，1) REM ψj1の解ﾍﾞｸﾄﾙzの中での位置
jjj1 = BN_Z(j1，2) REM qj1の解ﾍﾞｸﾄﾙzの中での位置
jj2 = BN_Z(j2，1) REM ψj2の解ﾍﾞｸﾄﾙzの中での位置
jjj2 = BN_Z(j2，2) REM qj2の解ﾍﾞｸﾄﾙzの中での位置
A(i, jj1) = A(i, jj1) - B1
A(i, jjj1) = A(i, jjj1) + H1
A(i, jj2) = A(i, jj2) - B2
A(i, jjj2) = A(i, jjj2) + H2
Next k
Next i
REM 領域要素番号kで領域積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
h = 0
For k = 1 to R_Count
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 特異点iと要素節点情報から積分値w0を計算
h = h + w0
Next k
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
z(i) = h
Next i
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 行番号から節点番号を特定する
j = i - BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
If BN_BP(j，1) = 1 Then
REM ψjが既知が場合
jj = BN_Z(j，1) REM ψjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM ψjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，1) REM ψjの値
ElseIf BN_BP(j，2) = 1 Then
REM qjが既知が場合
jj = BN_Z(j，2) REM qjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM qjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，2) REM qjの値
End If
Next i
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ﾍﾞｸﾄﾙzを求める
REM [出力部]
REM BN_Z(j，2) に従いz(j)から値を読み出し、出力
For j = 1 to BN_Count
jj = BN_Z(j，1) REM zの中でのψjの位置
REM z(jj)をψjの出力ﾙｰﾃｨﾝへ渡す
Next i
For j = 1 to BN_Count
jj = BN_Z(j，2) REM zの中でのqjの位置
REM z(jj)をqjの出力ﾙｰﾃｨﾝへ渡す
Next i
タグ：数値解析
