２次元電位分布の数値計算 の変更点
更新- 追加された行はこの色です。
- 削除された行はこの色です。
- サイエンス/２次元電位分布の数値計算 へ行く。
- サイエンス/２次元電位分布の数値計算 の差分を削除
[[公開メモ]]
#contents
* 概要 [#z820172d]
与えられた境界条件を満たす２次元電位分布を数値計算により求めたい。
* 基礎方程式 [#h7465ce4]
** 電位決定の基礎方程式 [#l6ddd6ff]
物質中におけるマクスウェル方程式において
&math(
\begin{cases}
\displaystyle\bm\nabla\times\bm E+\frac{\PD\bm B}{\PD t}=\bm o\\
\displaystyle\bm\nabla\cdot\bm B=0\\
\displaystyle\bm\nabla\times\bm H-\frac{\PD\bm D}{\PD t}=\bm i\\
\displaystyle\bm\nabla\cdot\bm D=\rho
\end{cases}
);
&math(\bm B=\bm 0); を仮定すれば、
&math(
\begin{cases}
\bm\nabla\times\bm E=\bm o\\
\displaystyle -\frac{\PD\bm D}{\PD t}=\bm i\\
\bm\nabla\cdot\bm D=\rho
\end{cases}
); (*)
第１式より電位の存在が導かれ、&math(\bm E=-\bm\nabla\phi);
&math(\bm D=\varepsilon\bm E); とともに第３式に代入すれば、
&math(
\bm\nabla\cdot\big(\varepsilon\bm\nabla\phi)=-\rho
);
特に、真電荷のない領域では
&math(
\bm\nabla\cdot\big(\varepsilon\bm\nabla\phi)=0
);
これが基礎方程式となる。
電場の z 成分がゼロのとき &math(E_z=0);、成分に分けて書けば
&math(\frac{\PD}{\PD x}\left(\varepsilon_x\frac{\PD}{\PD x}\phi\right)+
\frac{\PD}{\PD y}\left(\varepsilon_y\frac{\PD}{\PD y}\phi\right)=\rho);
媒質が等方的 (&math(\varepsilon_x=\varepsilon_y=\varepsilon);) かつ
均質 (&math(\varepsilon(\bm x)=\varepsilon);) であれば、
&math(\frac{\PD^2}{\PD x^2}\phi+
\frac{\PD^2}{\PD y^2}\phi=\triangle\phi=\rho/\varepsilon);
いわゆるポアソン方程式（&math(\rho=0); ならラプラス方程式）に帰着する。
** シート状伝導体における基礎方程式 [#o7161a92]
(*) 式の第２、第３式から電荷保存則
&math(-\bm\nabla\cdot\frac{\PD\bm D(\bm x,t)}{\PD t}=\bm\nabla\cdot\bm i=-\frac{d}{dt}\rho);
が導かれる。
定常状態では &math(\frac{d}{dt}\rho=0); より、
&math(\bm\nabla\cdot\bm i=0);
ここにシート状伝導体におけるオームの法則 &math(\bm i=\sigma\bm E=\begin{pmatrix}\sigma_x\\&\sigma_y\\&&0\end{pmatrix}\bm E);
電場と電位の関係式 &math(\bm E=-\bm\nabla\phi);
を代入すると、
&math(\bm\nabla\cdot(\sigma\bm\nabla\phi)=0);
これがシート状伝導体における電位の基礎方程式となる。
成分に分けて書けば、
&math(\frac{\PD}{\PD x}\left(\sigma_x\frac{\PD}{\PD x}\phi\right)+
\frac{\PD}{\PD y}\left(\sigma_x\frac{\PD}{\PD y}\phi\right)=0);
こちらも、コンダクタンスが等方的かつ一様であればポアソン方程式になる。
** 中心軸対称な電場分布に対する円柱座標の基礎方程式 [#s0e835b1]
基礎方程式
&math(\bm\nabla\cdot\big(\varepsilon\bm\nabla\phi)=-\rho);
を円柱座標
&math(
\begin{cases}
x=r\cos\theta\\
y=r\sin\theta\\
z=z\\
\end{cases}
);
で書き表すと、
&math(
\begin{cases}
\vspace{2mm}\displaystyle(\bm\nabla\phi)_r=\frac{\PD}{\PD r}\phi\\
\vspace{2mm}\displaystyle(\bm\nabla\phi)_\theta=\frac{1}{r}\frac{\PD}{\PD \theta}\phi\\
\displaystyle(\bm\nabla\phi)_z=\frac{\PD}{\PD z}\phi\\
\end{cases}
);
&math(
\bm\nabla\cdot\bm v=\frac{1}{r}\frac{\PD}{\PD r}(rv_r)+\frac{1}{r}\frac{\PD}{\PD\theta}v_\theta+\frac{\PD}{\PD z}v_z
);
より、
&math(
\bm\nabla\cdot(\varepsilon\bm\nabla\phi)
=\frac{1}{r}\frac{\PD}{\PD r}\left(r\varepsilon_x\frac{\PD}{\PD r}\right)\phi
+\frac{\PD}{\PD z}\left(\varepsilon_y\frac{\PD}{\PD z}\phi\right)
=-\rho
);
この場合には第1項の中に現れる r のため、単純なポアソン方程式にはならず、
方程式は座標 r に依存する
* ラプラス方程式を差分方程式に直す [#peb5f902]
２次元空間を &math(\Delta=\Delta x=\Delta y); を単位として格子状に分けて、
&math(
\phi_{p,q}=\phi(p\Delta x,q\Delta y)
);
とする。&math(\Delta); が十分に小さい時、
&math(\phi_{p+1,q}-\phi_{p,q}=\frac{\PD\phi}{\PD x}\Delta);
などが成り立つことを利用すると、ラプラス方程式を差分方程式に直すことができて、
&math(
\underbrace{\Big(\big(\phi_{p+1,q}-\phi_{p,q}\big)/\Delta-\big(\phi_{p,q}-\phi_{p-1,q}\big)/\Delta\Big)/\Delta}_{\displaystyle\PD^2\phi/\PD x^2}+
\underbrace{\Big(\big(\phi_{p,q+1}-\phi_{p,q}\big)/\Delta-\big(\phi_{p,q}-\phi_{p,q-1}\big)/\Delta\Big)/\Delta}_{\displaystyle\PD^2\phi/\PD y^2}=0
);
&math((\phi_{p+1,q}-2\phi_{p,q}+\phi_{p-1,q})+
(\phi_{p,q+1}-2\phi_{p,q}+\phi_{p,q-1})=0
);
&math(
\phi_{p,q}=\frac{1}{4}\big(\phi_{p+1,q}+\phi_{p-1,q}+\phi_{p,q+1}+\phi_{p,q-1}\big)
);
を得る。
つまり、この差分方程式を満たすような &math(\phi_{p,q}); は、元の方程式の近似的な解となる。
* ヤコビ法 [#r7e32476]
前項の最後に得た式を漸化式に置き換える。すなわち、
&math(
\phi_{p,q}^{(k+1)}=\frac{1}{4}\big(\phi_{p+1,q}^{(k)}+\phi_{p-1,q}^{(k)}+\phi_{p,q+1}^{(k)}+\phi_{p,q-1}^{(k)}\big)
);
何らかの方法により試行的な電位分布 &math(\phi_{p,q}^{(1)}); を定め、
上記の漸化式で分布を更新していった結果、もしこれが &math(\phi_{p,q}^{(\infty)}); に収束したならば、
&math(
\phi_{p,q}^{(\infty)}=\frac{1}{4}\big(\phi_{p+1,q}^{(\infty)}+\phi_{p-1,q}^{(\infty)}+\phi_{p,q+1}^{(\infty)}+\phi_{p,q-1}^{(\infty)}\big)
);
が成り立ち、この収束値が近似的に求めた電場分布を与えることになる。
このような手法はヤコビ法と呼ばれ、&math(\phi_{p,q}); の係数が比較的大きい場合にはよく収束することが知られている。
LANG:C
double calc_new_value(double t, double r, double b, double l)
{
return 0.25 * (t + r + b + l); // Laplace
}
参考：https://ja.wikipedia.org/wiki/%E3%83%A4%E3%82%B3%E3%83%93%E6%B3%95
* １つ飛ばしのヤコビ法 [#fde95cc5]
上記の漸化式では &math(\phi_{p,q}); はその上下左右の４つの格子点に依存しているが、
その他の点には依存していない。
すなわち、このような格子を考えれば、
&ref(checker-board.png);
赤いマスは白いマスに、白いマスは赤いマスにしか依存していない。
そこで、すべての格子点を一度に漸化式にかけるのではなく、
赤いマスと白いマスと交互に漸化式で更新してやれば、
&math(k); 回目と &math(k+1); 回目とに別々の記憶領域を取る必要がなくなる。
しかも恐らく、普通の２倍くらいの収束速度が得られそう？？？
* 境界条件 [#nda1dfaf]
領域の端では &math(\phi_{p+1,q}); などが領域外に出てしまうため、
何か条件を定めないと計算ができない。
領域の端を固定端とするには、領域外の値としてその固定値を入れてやれば良い。
領域の端を自由端とするには、領域外の値として &math(\phi_{p,q}); と同じ値を入れてやれば、
&math(\PD\phi/\PD x=0); などの条件を満たすことができる。
* 収束条件 [#zc07e772]
更新前と後とを比べて、変化が小さくなれば収束したと見なす。
具体的には、もっとも変化の大きかった格子点での変化量が規定値を下回れば収束したことにする。
もっといい方法があるかもしれないけれど、とりあえず、これで。
* 加速・減速 [#m1f004de]
通常のヤコビ法に比べて収束を早くしたければ、
LANG:c
new_value = calc_new_value();
new_value = old_value + (new_value - old_value) * speed;
のように、変化量を線形に増加 (1 < speed) させたり、減少 (0 < speed < 1) させたりすることも効果がある場合がある。
計算手法上は上記のままでは new_value - old_value の部分で桁落ちが心配されるため、
LANG:c
new_value = calc_new_value();
new_value = (1-speed) * old_value + speed * new_value;
と計算する方が良い。
このような方法を SOR 法と呼ぶ。https://ja.wikipedia.org/wiki/SOR%E6%B3%95
- 加速させると、うまく行けば収束が速くなるが、発散しやすくなる。
- 減速させると、通常のヤコビ法 (speed = 1) では発散してしまう場合にも解を得られる可能性がある。
* コードの概要 [#h48c7db6]
LANG:C
// n x n の正方領域を考える
// even_or_odd = 1 or 0 が与えられる
// pixel(p,q) で現在の値を得る
// set_pixel(p,q,v) で値を更新する
// mask(p,q) で固定端を判別する
// 更新差分の相対値の最大値を diff に残す
double half_step(int even_or_odd)
{
int p, q, nm1 = n-1;
int pstart = 1 + (even_or_odd ^ (q % 2));
double result = 0;
for(q = 1; q < nm1; q += 1 ){
for(p = pstart; p < nm1; p += 2 ){
if(mask(p, q)) continue; // 固定端は更新しない
double value = pixel(p, q);
if(isnan(value)) continue; // 自由端も更新しない
double new_value = calc_new_value(
replace_nan(pixel(p , q-1), value),
replace_nan(pixel(p+1, q ), value),
replace_nan(pixel(p , q+1), value),
replace_nan(pixel(p-1, q ), value)
);
// 速度を調整しつつ値を更新
set_pixel(p, q, (1-speed) * value + speed * new_value);
// 変化量を記録
double diff = fabs(new_value-value);
if( !isnan(diff) && diff > result ) result = diff;
}
}
return result;
}
double replace_nan(double value, double replace)
{
return isnan(value) ? replace : value;
}
* 徐々に粒度を上げる [#r7a997a4]
精確な電位を求めるには格子間隔 &math(\Delta); を小さくするのが望ましいけれど、
初めから小さくすると精確な解が得られるまでに多大な時間がかかってしまう。
そこで、細かい粒度で計算を始めるにあたり、粗い粒度で計算した結果を初期値として使うのが賢いやり方になる。
すなわち、まず粗い粒度で計算しておいて、徐々に粒度を上げていくのである。
** 例：一様な抵抗を持つ２次元伝導体への電流注入 [#i2bdfa72]
- 計算領域は ｘ, ｙ ともに±１の範囲 ＝ 一辺２の正方領域
- (±0.5, 0.0) を中心に、半径0.05の円形領域から電流を注入した
-- 電位を±１とした
*** 33pixel [#z5ad61f1]
&ref(sheet-33pixel.png);
*** 65pixel [#k71a7624]
&ref(sheet-65pixel.png);
*** 129pixel [#c3c7320b]
&ref(sheet-129pixel.png);
*** 257pixel [#u021c108]
&ref(sheet-257pixel.png);
* 計算ライブラリ [#cf949185]
計算効率はあまり追求しない方針で。
** pdesolver2d.h [#af873f8e]
LANG:cpp
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <float.h>
#include <algorithm>
#include <time.h>
#include <signal.h>
#ifndef __PDESOLVER2D
#define __PDESOLVER2D
// ２Ｄ正方領域で偏微分方程式 (PDE: Partial Differential Equation) を
// 解くためのクラス
//
// 領域サイズは n ピクセル x n ピクセル
// pixel に NaN を書き込むと、そこは自由端になる
// mask に非ゼロを書き込むと、そこは固定端になる
// 領域の周囲は必ず自由端もしくは固定端にしておかなければならない
//
// PDESolver2D 自体は抽象クラスなので、子クラスにて calc_new_value
// を実装して、周囲４点の値から新しい値を求める漸化式を与えないと
// 実体化できない
//
class PDESolver2D
{
protected:
int n;
double *pixels; // NaN を代入すると自由端になる
int *masks; // 非ゼロを代入すると固定端になる
int interrupted; // 強制終了する
double speed; // 加速、減速定数 (speed > 1 で加速)
public:
PDESolver2D(int n); // n x n の正方領域を確保する
virtual ~PDESolver2D();
int size() { return n; }
double pixel(int p, int q)
{ return pixels[n*q+p]; }
void pixel(int p, int q, double v)
{ pixels[n*q+p] = v; }
int mask(int p, int q)
{ return masks[n*q+p]; }
void mask(int p, int q, int v)
{ masks[n*q+p] = v; }
void clear() { memset(pixels, 0, n*n*sizeof(double));
memset(masks, 0, n*n*sizeof(int)); }
void interrupt(){ interrupted = 1; }
void set_speed(double s) { speed = s; }
// 全ピクセルで更新差分の相対値が delta 未満になるまで繰り返し計算する
// 更新差分の相対値とは更新差分をピクセル値で割った値のこと
// 計算が終了すれば非ゼロを返す
// 途中でキャンセルされればゼロを返す
// log がゼロでなければ stderr に表示する進行状況を log にも出力する
int solve(double delta, FILE *log = NULL);
// 指定されたファイルにデータを書き込む
void dump(const char *file_name);
// 現在の計算値を引き継いで２倍の解像度の領域を確保する
virtual PDESolver2D *create_double(int delete_this = 0);
protected: // 子クラスで実装すべき virtual 関数
// 自身の値 (c) と上下左右 (t, b, l, r) と座標値 (p,q) より、新しい値を計算して返す
virtual double calc_new_value(double c, double t, double r, double b, double l, int p, int q) = 0;
// サイズ n のインスタンスを作成して返す
virtual PDESolver2D *create_instance(int n) = 0;
protected: // 内部で用いるヘルパー関数
// create_double 内から呼ばれる。２倍領域のピクセル座標を
// 元画像に与えると、元画像のピクセル値を補完した値を返す
double interpolate(int p, int q);
// 自由端処理のためのヘルパー関数
double replace_nan(double value, double replace)
{ return isnan(value) ? replace : value; }
double replace_nan(int p, int q, double replace)
{ return replace_nan(pixel(p, q), replace); }
// 進行状況を表示する
void log_step(int i, double diff, time_t started_at, FILE *log);
// 奇数ピクセルあるいは偶数ピクセルだけを計算する
// ピクセル値の更新差分の相対値の最大値を返す
double half_step(int even_or_odd);
// 計算を１回分進める
// 内部では half_step を２回呼ぶ
// ピクセル値の更新差分の相対値の最大値を返す
double step();
};
#endif
** pdesolver2d.cpp [#a0e3d300]
LANG:cpp
#include "pdesolver2d.h"
PDESolver2D::PDESolver2D(int n)
{
speed = 1;
this->n = n;
interrupted = 0;
pixels = new double[n*n];
masks = new int[n*n];
clear();
}
PDESolver2D::~PDESolver2D()
{
delete pixels;
delete masks;
}
int PDESolver2D::solve(double delta, FILE *log)
{
time_t started_at = time(NULL);
int i;
double diff = -1;
int how_often = 53;
for(i = 1; !interrupted; i++) {
diff = step();
if(log && (i % how_often == 0))
log_step(i, diff, started_at, log);
if(diff < delta) break;
}
if(log) log_step(i, diff, started_at, log);
return !interrupted;
}
void PDESolver2D::log_step(int i, double diff, time_t started_at, FILE *log)
{
if(log!=stderr)
fprintf(log, "step (%6d, %e, %5.1fs)\n", i, diff, difftime(time(NULL), started_at));
fprintf(stderr, "step (%6d, %e, %5.1fs)\r", i, diff, difftime(time(NULL), started_at));
}
double PDESolver2D::half_step(int even_or_odd)
{
int p, q, nm1 = n-1;
double result = 0;
for(q = 1; q < nm1; q += 1 ){
for(p = 1 + (even_or_odd ^ (q % 2)); p < nm1; p += 2 ){
if(mask(p, q)) continue; // 固定端
double value = pixel(p, q);
if(isnan(value)) continue; // 自由端
double new_value = calc_new_value(
value,
replace_nan(p , q-1, value),
replace_nan(p+1, q , value),
replace_nan(p , q+1, value),
replace_nan(p-1, q , value),
p, q
);
pixel(p, q, value + speed * (new_value-value));
double diff = fabs(new_value-value)/(fabs(new_value)+fabs(value));
if( !isnan(diff) && diff > result ) result = diff;
}
}
return result;
}
double PDESolver2D::step()
{
return std::max( half_step(0), half_step(1) );
}
PDESolver2D *PDESolver2D::create_double(int delete_this)
{
int p, q;
PDESolver2D *result = create_instance(2*(n-1)+1);
for(p = 0; p < result->size(); p++){
for(q = 0; q < result->size(); q++){
result->pixel(p, q, interpolate(p, q));
}
}
if(delete_this) delete this;
return result;
}
double PDESolver2D::interpolate(int p, int q)
{
double i = 0;
double result = 0, pix;
int p2 = p/2, q2 = q/2;
if(p % 2){ // need interpolation for p
if(q % 2){ // need interpolation for q
if (!isnan(pix = pixel(p2 , q2 )) ) { result += pix; i++; }
if (!isnan(pix = pixel(p2+1, q2 )) ) { result += pix; i++; }
if (!isnan(pix = pixel(p2 , q2+1)) ) { result += pix; i++; }
if (!isnan(pix = pixel(p2+1, q2+1)) ) { result += pix; i++; }
return i > 0 ? result / i : NAN;
} else { // not need interpolation for q
if (!isnan(pix = pixel(p2 , q2 )) ) { result += pix; i++; }
if (!isnan(pix = pixel(p2+1, q2 )) ) { result += pix; i++; }
return i > 0 ? result / i : NAN;
}
} else { // not need interpolation for p
if(q % 2){ // need interpolation for q
if (!isnan(pix = pixel(p2 , q2 )) ) { result += pix; i++; }
if (!isnan(pix = pixel(p2 , q2+1)) ) { result += pix; i++; }
return i > 0 ? result / i : NAN;
} else { // not need interpolation for q
return pixel(p2, q2);
}
}
}
void PDESolver2D::dump(const char *file_name)
{
FILE *f = fopen(file_name, "w");
fwrite(pixels, sizeof(double), n*n, f);
fclose(f);
}
* 利用コード [#d3320029]
** main.cpp [#x0404645]
LANG:cpp
#include "pdesolver2d.h"
class SheetConductorSolver : public PDESolver2D
{
public:
SheetConductorSolver(int n)
: PDESolver2D(n)
{ }
PDESolver2D *create_instance(int n)
{ return new SheetConductorSolver(n); }
double calc_new_value(double c, double t, double r, double b, double l, int p, int q);
};
double SheetConductorSolver::calc_new_value(double c, double t, double r, double b, double l, int p, int q)
{
return 0.25 * (t + r + b + l); // Laplace
}
////////////////////////////////////////////////////////
double sqr(double v) { return v*v; }
void set_boundary(PDESolver2D *solver)
{
for(int p = 0; p < solver->size(); p++){
for(int q = 0; q < solver->size(); q++){
// outer boundary
if( p==0 || q==0 || p==solver->size()-1 || q==solver->size()-1 )
solver->pixel(p, q, NAN); // free edge
// current injection
double x = 2.0 * p / ( solver->size() - 1 ) - 1.0;
double y = 2.0 * q / ( solver->size() - 1 ) - 1.0;
// circular electrodes
if ( sqr(x+0.5) + sqr(y) < sqr(0.05) ) {
solver->mask(p, q, 1); // fixed edge
solver->pixel(p, q, -1.0);
}
if ( sqr(x-0.5) + sqr(y) < sqr(0.05) ) {
solver->mask(p, q, 1); // fixed edge
solver->pixel(p, q, +1.0);
}
}
}
}
PDESolver2D *solver = NULL;
// ^C が押されたら現在の計算を中断する
void signal_received(int sig)
{
if(solver)
solver->interrupt();
}
int main(int argc, const char *argv[])
{
int initial_size = 17;
solver = new SheetConductorSolver(initial_size);
// ^C が押されたら signal_received へ飛ぶ
signal(SIGINT, signal_received);
// 徐々に解像度を上げながら計算する
char fname[256];
double delta = 1e-6;
FILE *log;
// 0 1 2 3 4 5 6
// 17, 33, 65, 129, 257, 513, 1025
for(int i = 0; i < 7; i++) {
solver->set_speed(1);
set_boundary(solver);
sprintf(fname, "solver%dbefore.data", i);
solver->dump(fname);
sprintf(fname, "solver%d.log", i);
log = fopen(fname, "w");
int success = solver->solve(delta, log);
fclose(log);
if(success){
sprintf(fname, "solver%dafter.data", i);
} else {
sprintf(fname, "solver%dinterrupted.data", i);
}
solver->dump(fname);
fprintf(stderr, "\n*********** writing %s \n", fname);
int delete_original = 1;
solver = solver->create_double(delete_original);
}
}
** コンパイル＆実行 [#p1db3205]
LANG:console
$ ls
main.cpp pdesolver2d.cpp pdesolver2d.h
$ g++ -O4 main.cpp pdesolver2d.cpp
$ ls
*a.out main.cpp pdesolver2d.cpp pdesolver2d.h
$ ./a.out
step ( 337, 9.730739e-07, 0.0s)
*********** writing solver0after.data
step ( 948, 9.985805e-07, 0.0s)
*********** writing solver1after.data
step ( 3011, 9.996066e-07, 1.0s)
*********** writing solver2after.data
step ( 5274, 9.995926e-07, 1.0s)
*********** writing solver3after.data
step ( 464, 9.984991e-07, 1.0s)
*********** writing solver4after.data
step ( 958, 9.997507e-07, 5.0s)
*********** writing solver5after.data
step ( 1402, 9.997880e-07, 29.0s)
*********** writing solver6after.data
$
この系は素性がいいので、speed = 1.5 にするとさらに 1.5 倍速くなるけれど、
speed = 2 では発散する。
** 計算結果ファイル [#o224c208]
solver?after.data
のような名前のファイルが格段の計算結果となる。
LANG:cpp
double pixels[n*n];
の内容を単にダンプしたものなので、sqrt( (ファイルサイズ)/8バイト ) で n を求めて、
n*n の二次元配列として double 値を読み込み、グラフソフトで表示することになる。
Igor Pro なら、
LANG:igor
function show(filename)
string filename;
GBLoadWave/B/T={4,4}/W=1/O filename
variable n = sqrt(numpnts(wave0))
make/O/N=(n,n) image
wave wave0, image
image = wave0[p+q*n]
setscale x,0,2,image
setscale y,-1,1,image
end
を使って、
LANG:command
*show("z:/electric-field2/solver6after.data")
*newimage image
*ModifyImage image ctab= {*,*,Rainbow,0}
とすればカラー表示できる。
&ref(sheet-conductor1.png);
** 等高線 [#zeac2a80]
[グラフ]-[グラフに追加]-[等高線プロット]
から、
[アピアランス] の
- [ラインカラー] を [すべて同じカラー]-[白]
- [ラベル詳細] で、
-- [背景] を [透明]
-- [ラベルテキスト]-[カラー]-[黒]
とすれば、
LANG:console
*AppendMatrixContour/T image
*ModifyContour image rgbLines=(65535,65535,65535),labelBkg=1,labelRGB=(0,0,0)
と同等の動作になり、
&ref(sheet-conductor2.png);
が得られる。
** 電場強度を見るには [#wdc84a72]
強度のみであれば &math((\phi_{p+1,q}-\phi_{p,q})^2+(\phi_{p,q+1}-\phi_{p,q})^2); をプロットすればいい。
Igor Pro なら、
LANG:console
*duplicate image, image2
*image2 = sqrt((p == 0 ? image[p+1][q]-image[p][q] : image[p][q]-image[p-1][q])^2+(q == 0 ? image[p][q+1]-image[p][q] : image[p][q]-image[p][q-1])^2)
*image2 = image2 == 0 ? NaN : image2
** 電場の方向を見るには [#t3ec705f]
function calcfield()
wave image
variable n = dimsize(image,0)
variable d = 20
variable pp, qq, nn = trunc((n-1)/d)
variable ma = 0
make/O/N=(nn*nn) posx, posy
make/O/N=(nn*nn,2) field
for(pp=0; pp<nn-1; pp+=1)
for(qq=0; qq<nn-1; qq+=1)
posx[pp+qq*nn] = DimOffset(image, 0) + DimDelta(image, 0) * pp * d
posy[pp+qq*nn] = DimOffset(image, 1) + DimDelta(image, 1) * qq * d
variable fieldx = -image[pp*d+1][qq*d ]+image[pp*d][qq*d]
variable fieldy = image[pp*d ][qq*d+1]-image[pp*d][qq*d]
variable fieldr = sqrt(fieldx^2+fieldy^2)
if (fieldr > ma)
ma = fieldr
endif
field[pp+qq*nn][1] = atan2(fieldy, fieldx)
endfor
endfor
field[][0] = 40 // 長さ
end
として、
LANG:console
*calcfield()
*appendtograph/t/l posy vs posx
*ModifyGraph mode=3,arrowMarker(posy)={field,0.5,0,0.5,1},msize=0.1
*ModifyGraph rgb(posy)=(0,0,0)
などとすれば、下記のようなグラフを得る。
&ref(sheet-conductor3.png);
* 円柱座標の場合の計算 [#jcaaa019]
誘電率が均一な場合、上記のラプラシアンを変形して、
&math(
\Delta\phi
&=\frac{1}{r}\frac{\PD}{\PD r}\left(r\frac{\PD}{\PD r}\right)\phi
+\frac{\PD^2}{\PD z^2}\phi\\
&=\frac{1}{r}\frac{\PD}{\PD r}\phi
+\frac{\PD^2}{\PD r^2}\phi
+\frac{\PD^2}{\PD z^2}\phi\\
&=0
);
そのまま差分方程式に直すと、
&math(
(\phi_{p+1,q}-\phi_{p-1,q})/2p\Delta r^2
+(\phi_{p+1,q}-2\phi_{p,q}+\phi_{p-1,q})/\Delta r^2
+(\phi_{p,q+1}-2\phi_{p,q}+\phi_{p,q-1})/\Delta r^2=0
);
整理すると、
&math(
\phi_{p,q}=(\phi_{p+1,q}-\phi_{p-1,q})/8p
+(\phi_{p+1,q}+\phi_{p-1,q}+\phi_{p,q+1}+\phi_{p,q-1})/4
);
calc_new_value 中でこの式を使えば正しい計算結果が得られる。
すなわち、
LANG:cpp
class CylindricalSolver : public PDESolver2D
{
public:
CylindricalSolver(int n)
: PDESolver2D(n)
{ }
PDESolver2D *create_instance(int n)
{ return new CylindricalSolver(n); }
double calc_new_value(double c, double t, double r, double b, double l, int p, int q);
};
double CylindricalSolver::calc_new_value(double c, double t, double r, double b, double l, int p, int q)
{
return 0.25 * (t + r + b + l) + 0.125 * (r - l) / p; // Cylindrical Laplace
}
とすればソルバーの完成。
あとは、
LANG:cpp
void set_boundary(PDESolver2D *solver, double r1, double r2, double d)
{
for(int p = 0; p < solver->size(); p++){
for(int q = 0; q < solver->size(); q++){
// outer boundary
if( p==0 || q==0 || p==solver->size()-1 || q==solver->size()-1 )
solver->pixel(p, q, NAN); // free edge
// coordinate
double r = 2.0 * p / ( solver->size() - 1 );
double z = 2.0 * q / ( solver->size() - 1 ) - 1.0;
if(d/2<z && z<d/2+0.1 && r<r1) {
solver->pixel(p, q, 0); // ground
}
if(-d/2<z && z<-d/2-0.1 && r<r2) {
solver->pixel(p, q, 1); // biased
}
}
}
}
のようにして境界条件を与えればいい。
* 誘電率が不連続に変化する場合 - 有限要素法を使う [#ye06460c]
上記のように誘電率や伝導度が均一の時は差分方程式をそのままヤコビ法で解いてうまく行くのに対して、
誘電率が不連続に変化する系に対して無理やり差分方程式を立てて解こうとしても、
繰り返し計算が発散してしまいうまく計算できない。
参考：~
http://www.yamamo10.jp/yamamoto/study/electromagnetics/change_electric_constant/laplace.pdf
> 差分法はこのように誘電率が異なる静電場の計算には適さない
そこでどうするかというと、有限要素法により漸化式を得るのが良いそうだ。
** 平面での漸化式 [#xab04648]
停留条件は、
&math(
\frac{\PD U}{\PD \phi_{p,q}}
&=\frac{\varepsilon_{p-1/2,q-1/2}}{2}\big(2\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q-1}\big)\\
&+\frac{\varepsilon_{p+1/2,q-1/2}}{2}\big(2\phi_{p,q}-\phi_{p+1,q}-\phi_{p,q-1}\big)\\
&+\frac{\varepsilon_{p-1/2,q+1/2}}{2}\big(2\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q+1}\big)\\
&+\frac{\varepsilon_{p+1/2,q+1/2}}{2}\big(2\phi_{p,q}-\phi_{p+1,q}-\phi_{p,q+1}\big)\\
&=0
);
これを整理した漸化式は、
&math(
&\big[2(\varepsilon_{p-1/2,q-1/2}+\varepsilon_{p+1/2,q-1/2}+\varepsilon_{p-1/2,q+1/2}+\varepsilon_{p+1/2,q+1/2})\big]\phi_{p,q}\\
&=\varepsilon_{p-1/2,q-1/2}\big(\phi_{p-1,q}+\phi_{p,q-1}\big)+\varepsilon_{p+1/2,q-1/2}\big(\phi_{p+1,q}+\phi_{p,q-1}\big)\\
&\ +\varepsilon_{p-1/2,q+1/2}\big(\phi_{p-1,q}+\phi_{p,q+1}\big)+\varepsilon_{p+1/2,q+1/2}\big(\phi_{p+1,q}+\phi_{p,q+1}\big)
);
となる。
これも上下左右の格子点の値しか使っていないため、加えて epsilon(x,y) を与えてやることで
上記のソルバーを使って解ける。
** 円柱座標での漸化式 [#c0411a18]
円柱座標の場合、対応する停留条件は、
&math(
\frac{\PD U[\phi]}{\PD \phi_{p,q}}&=\frac{\pi \Delta}{2}\Big\{
p\varepsilon_{p+1/2,q+1/2}\big(-\phi_{p+1,q}+\phi_{p,q}-\phi_{p,q+1}+\phi_{p,q}\big)+\\
&\hspace{11mm}\ (p-1)\varepsilon_{p-1/2,q+1/2}\big(\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q+1}+\phi_{p,q}\big)+\\
&\hspace{11mm}\ p\varepsilon_{p+1/2,q-1/2}\big(-\phi_{p+1,q}+\phi_{p,q}+\phi_{p,q}-\phi_{p,q-1}\big)+\\
&\hspace{11mm}\ (p-1)\varepsilon_{p-1/2,q-1/2}\big(\phi_{p,q}-\phi_{p-1,q}+\phi_{p,q}-\phi_{p,q-1}\big)\Big\}\\
&=0
);
整理すると、
&math(
&2\big[p\varepsilon_{p+1/2,q+1/2}+(p-1)\varepsilon_{p-1/2,q+1/2}+p\varepsilon_{p+1/2,q-1/2}+(p-1)\varepsilon_{p-1/2,q-1/2}\big]\phi_{p,q}\\
&=p\varepsilon_{p+1/2,q+1/2}\big(\phi_{p+1,q}+\phi_{p,q+1}\big)+
(p-1)\varepsilon_{p-1/2,q+1/2}\big(\phi_{p-1,q}+\phi_{p,q+1}\big)+\\
&\hspace{5mm}p\varepsilon_{p+1/2,q-1/2}\big(\phi_{p+1,q}+\phi_{p,q-1}\big)+
(p-1)\varepsilon_{p-1/2,q-1/2}\big(\phi_{p-1,q}+\phi_{p,q-1}\big)\\
);
となる。
* コメント・質問 [#yf7ffa8e]
#article_kcaptcha
**自由端の境界条件についての質問 [#rb74db86]
>[[くろぐろ]] (&timetag(2022-11-12T09:29:59+09:00, 2022-11-12 (土) 18:29:59);)~
~
メモを拝見させていただきました. 等電位線のシミュレーションを行おうとしていたのでとても役立ちました. 表題の件について質問です.~
自由端において, 領域外の値としてφpqを使用することで∂φ/∂x=0が満たされるといのことですが, この条件は領域端を走る電気力線に対して等電位線が垂直に交わることを表しているのでしょうか? また, φpqを使用すると, なぜこの条件が満たされるのでしょうか? ~
この質問に答えていただけると幸いです.~
//
- 前提として、このように離散化した数値計算においては隣のセル同士の間で電位を差し引きしたものが電場になります。 ですので、領域外の電位を境界上の電位と等しくすれば、境界上のセルにおいて「隣のセルとの差分がゼロ」になり、これは面直方向の電場がゼロという条件に相当するのだと考えています。 -- [[武内（管理人）]] &new{2022-11-12 (土) 21:10:13};
- 返答ありがとうございます. 納得できました. -- [[くろぐろ]] &new{2022-11-13 (日) 18:55:13};
#comment_kcaptcha
Counter: 18402 (from 2010/06/03),
today: 1,
yesterday: 5