もっと高速な方法はあるでしょうか？

　　もっと高速な方法はあるでしょうか？ 島村1243 2007/12/13 14:21:55  (修正1回)
　　└続き 島村1243 2007/12/13 14:23:05  (修正1回)
　　　├(1)について（処理手順の改良） 山中和義 2007/12/14 08:55:07 
　　　│├!もったいない、整理すれば、４倍速。 SECOND 2007/12/14 21:34:04 
　　　││└確かに速いです。 島村1243 2007/12/14 22:33:52 
　　　││　└どっちでしょう？ SECOND 2007/12/14 22:50:11  (修正1回)
　　　│└２０秒速くなりました。 島村1243 2007/12/14 22:06:46 
　　　└(3)について各点での流線 山中和義 2007/12/15 09:43:05 
　　　　├!等ポテンシャルを色の濃淡で表現する 山中和義 2007/12/15 11:03:55 
　　　　└これを利用すると流れの方向がわかるので、... 山中和義 2007/12/15 20:11:13 
　　　　　└とても数理的で良いですね。 島村1243 2007/12/15 23:15:50  (修正1回)
　　　　　　└接地導体球の作画に問題が有る様 島村1243 2007/12/16 07:37:08 
　　　　　　　└対策案できました。 島村1243 2007/12/16 10:44:04  (修正1回)
　　　　　　　　└ポテンシャル関数f(x,y)なら、電位分布は等... 山中和義 2007/12/16 13:56:21 
　　　　　　　　　└完璧に目的達成です！ 島村1243 2007/12/16 19:04:03  (修正1回)
　　　　　　　　　　└プログラムの続き 島村1243 2007/12/16 19:05:45

　　もっと高速な方法はあるでしょうか？ 島村1243 2007/12/13 14:21:55  (修正1回) ツリーへ

もっと高速な方法はあるでしょうか？ 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/13 14:21:55 ** この記事は1回修正されてます

電気影像法（映像法）の理屈と等角写像を組み合わせて、複数の点電荷が有る場合の２次元静電場の画像を描か
せる方法を模索中です。

z平面の（a1,j0）と(a2,j0)の２点にそれぞれq1、q2の電荷（+,-を考慮）を置いた場合の等電位分布と電気
力線分布図は、
w=q1*Log(z-a1)+q2*Log(z-a2)、w=u+j*v （ただし左記以外の定数は省略）
をz=f(w)と解析的に表す事が出来れば電位分布は

j=sqr(-1)
For u=-10 to 10 step 0.2
For v=-pi to pi step pi/20
w=u+j*v
z=f(w)
plot lines:Re(z),Im(z);
next
plot lines
next

の様にして高速に、しかも綺麗に作図出来るはずです（電気力線分布も同様）。

ちなみにa2=-a1,q1=q2=+1の場合（２個の正電荷がy軸対称に設置された場合）は
w=Log(z-a1)+Log(z+a1) から z=f(w)=sqr(a1^2+exp(w))を得て作図が高速に行えました。

しかしa1<>a2,q1<>q2の様な一般的な場合は解析解が得られない（と言うより解析方法が分からない）ので
上記写像の理屈で作画が出来ません。

電位分布だけなら各点電荷による電位を合成した式に、以前の投稿「根軌跡図の作成」で山中様に教えて頂い
た作画方法
http://freebbs.around.ne.jp/article/b/basic/96/sfpfow/vzgpmu.html#vzgpmu
を利用して下記の様に作成できます。

　　└続き 島村1243 2007/12/13 14:23:05  (修正1回) ツリーへ

Re: もっと高速な方法はあるでしょうか？ 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/13 14:23:05 ** この記事は1回修正されてます

続き

!接地導体球と1個の点電荷による静電場の電位分布作図
!等式に基づくドット画像を描く山中氏の方法を利用

OPTION ARITHMETIC complex
LET Xa=-6
LET Xb=-Xa
LET Ya=-6
LET Yb=-Ya

SET WINDOW Xa,Xb,Ya,Yb !描画範囲
DRAW grid
ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算
LET c_Div=3 !分解数
LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る
LET dy=(Yb-Ya)/b/c_Div
LET cEps=1E-3 !誤差　※精度が気になる!?
LET j=SQR(-1) !虚数単位

LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標セット
LET a2=r^2/a1 !影像電荷のx座標計算
LET q1=1 !点電荷の大きさセット
LET q2=-r/a1*q1 !影像電荷の大きさ計算

!点電荷の位置に赤丸印を表示
SET AREA COLOR "red"
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

!接地導体球の位置に緑丸印を表示
SET AREA COLOR "green"
DRAW disk WITH SCALE(r)*SHIFT(0,0)

SET POINT STYLE 1
FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定
IF wb<=0 THEN
SET POINT COLOR "black"
ELSE
SET POINT COLOR "red"
END IF
FOR k=0 TO c_Div*b-1 !Ｙ座標　※ドットをもとに
LET y=Ya+dy*k ![Ya,Yb]
FOR i=0 TO c_Div*a-1 !Ｘ座標
LET x=Xa+dx*i ![Xa,Xb]
LET z=x+y*j !空間上の点を計算する
IF z<> a1 AND z<>a2 THEN
LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位
IF ABS(w-wb) < cEps THEN PLOT POINTS: x,y !点画
END IF
NEXT i
NEXT k
NEXT wb
END

しかし、
（１）1.7GHzのPCでは作図に長い時間を要する。
（２）ドット表示なので線が粗い。また、誤差1/1000設定でもドット幅が広い部分が生じてしまう。
（３）電気力線分布が得られない。
と言う難点があります。

何か良い解決策は有るでしょうか？

　　　├(1)について（処理手順の改良） 山中和義 2007/12/14 08:55:07  ツリーへ

Re: 続き 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/14 08:55:07

(1)について（処理手順の改良）

等高線の数だけ複素平面を走査して毎回wを算出すると時間を費やします。
FOR文のネストをかえると多少速くなります。

「１回の走査で該当する等高線を描く」ように改修する場合（FOR文の部分抜粋）

SET POINT STYLE 1
FOR k=0 TO c_Div*b-1 !Ｙ座標　※ドットをもとに
LET y=Ya+dy*k ![Ya,Yb]
FOR i=0 TO c_Div*a-1 !Ｘ座標
LET x=Xa+dx*i ![Xa,Xb]

LET z=x+y*j !空間上の点を計算する
IF z<>a1 AND z<>a2 THEN
LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位

FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定
IF wb<=0 THEN
SET POINT COLOR "black"
ELSE
SET POINT COLOR "red"
END IF
IF ABS(w-wb) < cEps THEN PLOT POINTS: x,y !点画
NEXT wb
END IF
NEXT i
NEXT k

　　　│├!もったいない、整理すれば、４倍速。 SECOND 2007/12/14 21:34:04  ツリーへ

Re: (1)について（処理手順の改良） 返事を書くノートメニュー

SECOND <jjqdmekgpt> 2007/12/14 21:34:04

!もったいない、整理すれば、４倍速。

OPTION ARITHMETIC complex
LET Xa=-6
LET Xb=-Xa
LET Ya=-6
LET Yb=-Ya

SET WINDOW Xa,Xb,Ya,Yb !描画範囲
DRAW grid
ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算
LET c_Div=3 !分解数
LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る
LET dy=(Yb-Ya)/b/c_Div
LET cEps=1E-3 !誤差　※精度が気になる!?
! LET j=SQR(-1) !虚数単位　　…使わない

LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標セット
LET a2=r^2/a1 !影像電荷のx座標計算
LET q1=1 !点電荷の大きさセット
LET q2=-r/a1*q1 !影像電荷の大きさ計算

!点電荷の位置に赤丸印を表示
SET AREA COLOR "red"
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

!接地導体球の位置に緑丸印を表示
SET AREA COLOR "green"
DRAW disk WITH SCALE(r)*SHIFT(0,0)

SET POINT STYLE 1
FOR y=Ya TO Yb STEP dy ![Ya,Yb]
FOR x=Xa TO Xb STEP dx ![Xa,Xb]
LET z=complex(x,y) !空間上の点z

IF z<>a1 AND z<>a2 THEN
LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位
FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定
IF ABS(w-wb) < cEps THEN
IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red"
PLOT POINTS: x,y !点画
END IF
NEXT wb
END IF

NEXT x
NEXT y

END

　　　││└確かに速いです。 島村1243 2007/12/14 22:33:52  ツリーへ

Re: !もったいない、整理すれば、４倍速。 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/14 22:33:52

確かに速いです。

SECONDさん、いつもご教示有難う御座います。
確かに速く、ストレスが感じなくなりますね。

SECONDさんのテクニックで４倍速になる理由が何なのかを知るために、山中さんのコードとの違いを
見たのですが

（１）z=u+j*vを使用せず、z=complex(x,y)を使用している。
（２）画面走査に「FOR k=　　」のループを使用せず、刻み幅設定は「STEP」で指定し、FORループを
２段に削減している。

でした。今回の場合、上記（１）と（２）のどちらが大きな効力を発揮しているのでしょうか？

　　　││　└どっちでしょう？ SECOND 2007/12/14 22:50:11  (修正1回) ツリーへ

Re: 確かに速いです。 返事を書くノートメニュー

SECOND <jjqdmekgpt> 2007/12/14 22:50:11 ** この記事は1回修正されてます

どっちでしょう？
実行頻度の最も多い　最も内側のステートメントの処理を軽くする
努力が、一番効果します。文の順番も重要です。
それから、
使用する変数の種類を減らす工夫をすると、不思議な事に、
自然に文の方が誘導してくれる傾向があります。

これは、山中氏のアイデアがよかったから、マグレですね。(^^;)

それより、逆写像の式を求める方が、大事ですね。私にもわかりません。

　　　│└２０秒速くなりました。 島村1243 2007/12/14 22:06:46  ツリーへ

Re: (1)について（処理手順の改良） 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/14 22:06:46

２０秒速くなりました。

山中さん、ご教示有難う御座います。
同じ全画面走査を何回も繰り返すよりも、全画面走査は1回の方が速くなると言う論理ですね。
私のPCでは原作が1分５５秒かかって完了に対し、山中さんの改修コード適用の場合は1分３５秒でした。

　　　└(3)について各点での流線 山中和義 2007/12/15 09:43:05  ツリーへ

Re: 続き 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/15 09:43:05

(3)について　各点での流線

!複素ポテンシャル

OPTION ARITHMETIC complex

LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ

DEF f(z)=q1*LOG(z-a1)+q2*LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み）

SET WINDOW -6,6,-6,6 !描画範囲
DRAW grid

SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

LET h=0.01 !増分

FOR x=-6 TO 6 STEP 0.5 !複素平面を走査する
FOR y=-6 TO 6 STEP 0.5
WHEN EXCEPTION IN
LET w=f(complex(x,y)) !ベクトル成分を得る　※流線
LET wx=f(complex(x+h,y)) !偏微分係数から
LET wy=f(complex(x,y+h))
LET vy=-Im(wx-w)/h !∂ψ/∂x
LET vx= Im(wy-w)/h !∂ψ/∂y

LET ww=complex(vx,vy) !成分に応じたベクトルを得る
DRAW arrow WITH ROTATE(arg(ww))*SCALE(ABS(ww)*0.7)*SHIFT(x,y) !※0.7 調整要
USE
END WHEN
NEXT y
NEXT x

PICTURE arrow !矢印を描く
PLOT LINES: 0,0; 1,0
PLOT LINES: 1,0; 0.8, 0.1
PLOT LINES: 1,0; 0.8,-0.1
END PICTURE

END

　　　　├!等ポテンシャルを色の濃淡で表現する 山中和義 2007/12/15 11:03:55  ツリーへ

Re: (3)について各点での流線 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/15 11:03:55

!等ポテンシャルを色の濃淡で表現する

OPTION ARITHMETIC complex

LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ

DEF f(z)=q1*LOG(z-a1)+q2*LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み）

SET WINDOW -6,6,-6,6 !描画範囲

SET POINT STYLE 1

FOR x=-6 TO 6 STEP 0.02 !複素平面を走査する
FOR y=-6 TO 6 STEP 0.02
WHEN EXCEPTION IN
LET w=f(complex(x,y))
!LET c=MOD(Re(w)*1000,256)/256 !電位　※1000 調整要
LET c=MOD(Im(w)*1000,256)/256 !電気力線
SET COLOR MIX(1) 0,0,c !青の濃淡
PLOT POINTS: x,y
USE
END WHEN
NEXT y
NEXT x

DRAW grid !座標

SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a2,0)

END

　　　　└これを利用すると流れの方向がわかるので、... 山中和義 2007/12/15 20:11:13  ツリーへ

Re: (3)について各点での流線 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/15 20:11:13

これを利用すると流れの方向がわかるので、タートル（亀）をそれに沿って移動させると、、、
（川の流れに身をまかせる、、、感じです。）

!タートルグラフィックスで電気力線を描く

OPTION ARITHMETIC complex

LET r=1 !導体球の半径
LET a1=4 !点電荷のx座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ

DEF f(z)=q1*LOG(z-a1)+q2*LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み）

SET WINDOW -6,6,-6,6 !描画範囲
DRAW grid

LET L=0.2 !移動量

LET h=0.001 !増分

FOR th=0 TO 360 STEP 5 !点電荷から放射状に
LET theta=RAD(th)
LET x=a1
LET y=0
PLOT LINES: x,y;

DO
WHEN EXCEPTION IN
LET w=f(complex(x,y)) !ベクトル成分を得る　※流線
LET wx=f(complex(x+h,y)) !偏微分係数から
LET wy=f(complex(x,y+h))
LET vy=-Im(wx-w)/h !∂ψ/∂x
LET vx= Im(wy-w)/h !∂ψ/∂y

LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る
USE
END WHEN

LET x=x+L*COS(theta) !移動させる
LET y=y+L*SIN(theta)
PLOT LINES: x,y;

IF ABS(x)>6 THEN EXIT DO !表示範囲外なら
IF ABS(y)>6 THEN EXIT DO

IF ABS(x*x+y*y)<r THEN EXIT DO !接地導体球内なら
LOOP

PLOT LINES !次へ
NEXT th

SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

END

　　　　　└とても数理的で良いですね。 島村1243 2007/12/15 23:15:50  (修正1回) ツリーへ

Re: これを利用すると流れの方向がわかるので、... 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/15 23:15:50 ** この記事は1回修正されてます

とても数理的で良いですね。

山中さん、いろいろなご検討有り難うございます。
電位分布（または電気力線）を色の濃淡に変換する案は、何故、連続的に変化している電位に縞模様の濃淡
が出きるのか？と言う点が理解できない状態ですが、視覚的にはとても素晴らしい発想ですね。プログラム
コードを良く見て何故？を考えようと思います。

最後に示して頂いたタートルグラフィックスは、数理的な扱いも分かり易く、かつ、当初求めていた描画の
目的が達成出来る解決策だと思います。山中さんのプログラム例では、正電荷から負電荷に向かって電気力
線が放出されると言う考え方で、１ヶ所からの電気力線放出のみを計算していました。
正電荷が２個有る場合にも、同じ様な考え方で２ヶ所について放出計算をすれば正解が得られるかを試行し
てみます。

さて、解決策を示して頂き胸わくわくだったのですが、パソコンの描画結果をよくよく眺めてみると、何か
おかしいなと気づきました。

１個の正点電荷とその近くに接地導体球が有る二次元空間の電位分布は、電気影像法に依れば正電荷とそれ
の影像電荷のそれぞれの電荷によって生じる電位分布を代数加算して求められます。
私が最初に提起したプログラムはその方法を使って電位分布を計算していますが、実際に計算される接地導
体球(半径１としました）の表面電位はゼロになり、かつ、接地導体球の周囲には球（完全球ではない）状の
電位分布が描かれているので正しい結果を表している事が分かります。

ところが、複素ポテンシャル関数w=q1*Log(z-a1)+q2*Log(z-a2)を使って山中さんが作成してくださった
電位分布状況図について接地導体球（原点を中心とした半径１）の周辺について見ると、等電位線が接地導
体球表面と交わる形に描かれています。これは大きさの有る接地導体球にかんしては正しい作図を表してい
ない（点電荷の電位分布図としては正しい作図の様に思えますが）と言う事です。

つまり、「接地導体球内に仮想された影像電荷と真の点電荷の全体の複素ポテンシャル関数をそれぞれの複素
ポテンシャル関数を単純に加算（例えばq1*Log(z-a1)+q2*Log(z-a2)の様に）して求め、その式から電位
分布・電気力線分布を作図すれば良い」と言う私の考え方は根本的に正しく無かったのではないか？と言うこ
とです。

　　　　　　└接地導体球の作画に問題が有る様 島村1243 2007/12/16 07:37:08  ツリーへ

Re: とても数理的で良いですね。 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/16 07:37:08

接地導体球の作画に問題が有る様

複素ポテンシャル関数はLog関数、一方原点を中心に半径１の接地導体円を描いた座標系はLogとは無関係。
その辺に誤りの根本が有るように感じています。

　　　　　　　└対策案できました。 島村1243 2007/12/16 10:44:04  (修正1回) ツリーへ

Re: 接地導体球の作画に問題が有る様 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/16 10:44:04 ** この記事は1回修正されてます

対策案できました。

山中さんがご教示くださったタートルグラフィックスの手法と、SECONDさんがご教示くださった高速
描画手法を組み合わせて、（複素ポテンシャル関数は諦めて）普通のポテンシャル関数を使って描画
したら、（電位分布は誤差の為に線に幅が生じていますが）正しい分布状態の作図に成功しましたの
で以下に示します。
お二方には感謝申し上げます。

!----点電荷と接地導体球の電場（電位分布と電気力線分布）作図----
!-------山中和義氏とSECOND氏の方法を合体----

LET Xa=-10
LET Xb=-Xa+4
LET Ya=Xa-2
LET Yb=-Ya+2
SET WINDOW Xa,Xb,Ya,Yb !描画範囲
DRAW grid
ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算
LET c_Div=3 !分解数
LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る
LET dy=(Yb-Ya)/b/c_Div
LET cEps=1E-3 !誤差　※精度が気になる!?

LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ
DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数

SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)
SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

!--電位分布の作図---
SET POINT STYLE 1
FOR y=Ya TO Yb STEP dy ![Ya,Yb]
FOR x=Xa TO Xb STEP dx ![Xa,Xb]
IF x<>a1 AND x<>a2 THEN
LET w=f(x,y) !点(x,y)における合成電位
FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定
IF ABS(w-wb) < cEps THEN
IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red"
PLOT POINTS: x,y !点画
END IF
NEXT wb
END IF
NEXT x
NEXT y

!----電気力線の描画----
LET L=0.2 !移動量
LET h=0.001 !増分
FOR th=0 TO 360 STEP 3 !点電荷から放射状に
LET theta=RAD(th)
LET x=a1
LET y=0
PLOT LINES: x,y;

DO
WHEN EXCEPTION IN
LET w=f(x,y) !電位を得る
LET wx=f(x+h,y) !偏微分係数から
LET wy=f(x,y+h)
LET vx=-(wx-w)/h !∂f/∂x
LET vy=-(wy-w)/h !∂f/∂y
LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る
USE
END WHEN

LET x=x+L*COS(theta) !移動させる
LET y=y+L*SIN(theta)
PLOT LINES: x,y;
IF ABS(x)>Xb THEN EXIT DO !表示範囲外なら
IF ABS(y)>Yb THEN EXIT DO
IF ABS(x*x+y*y)<r THEN EXIT DO !接地導体球内なら
LOOP
PLOT LINES !次へ
NEXT th

!---画面仕上げ用（点電荷位置・接地導体球）再描画----
SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)
SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

END

　　　　　　　　└ポテンシャル関数f(x,y)なら、電位分布は等... 山中和義 2007/12/16 13:56:21  ツリーへ

Re: 対策案できました。 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/16 13:56:21

ポテンシャル関数f(x,y)なら、電位分布は等高線アルゴリズムで描けます。

!関数f(x,y)の等高線を描く

LET cEps=0.001 !許容誤差
LET h=0.001 !増分（偏微分係数）

SUB newton(fc,x,y, fx,fy,grad2) !ニュートン法
DO
LET ff=f(x,y) !∇f
LET fx=(f(x+h,y)-ff)/h
LET fy=(f(x,y+h)-ff)/h
LET grad2=fx*fx+fy*fy

IF grad2<1e-10 THEN
LET x=1e30
EXIT SUB
END IF

LET t=(fc-ff)/grad2
LET x=x+t*fx
LET y=y+t*fy
LOOP WHILE t*t*grad2>cEps*cEps
END SUB

SUB contour(fc,x,y,d) !等高線を描く
LET i=0
DO
CALL newton(fc,x,y, fx,fy,grad2)
IF ABS(x)+ABS(y)>1e10 THEN EXIT SUB
IF i=0 THEN
PLOT LINES: x,y; !始点
LET x0=x
LET y0=y
ELSE
PLOT LINES: x,y; !折れ線でつなげる
END IF
IF i>2 AND (x-x0)^2+(y-y0)^2<d*d THEN EXIT DO !始点近傍なら、終了

LET t=d/SQR(grad2)
LET x=x+fy*t
LET y=y-fx*t

LET i=i+1
LOOP
PLOT LINES: x0,y0 !閉じる
END SUB
!---------- ここまでがサブルーチン

!SET WINDOW -6,6,-6,6 !表示領域
SET WINDOW -15,15,-15,15 !表示領域
DRAW grid !座標

LET r=1 !導体球の半径
LET a1=4 !点電荷のx座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ

DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数

SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)

SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

LET d=0.015 !増分
FOR fc=0 TO 0.5 STEP d !等高線の関数値
CALL contour(fc,1,1,d)
NEXT fc

FOR fc=0 TO 0.06 STEP d !※左端
CALL contour(fc,-10,1,d)
NEXT fc

END

　　　　　　　　　└完璧に目的達成です！ 島村1243 2007/12/16 19:04:03  (修正1回) ツリーへ

Re: ポテンシャル関数f(x,y)なら、電位分布は等... 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/16 19:04:03 ** この記事は1回修正されてます

完璧に目的達成です！

山中さん、幾度ものご教示有難う御座いました。
当初は複素ポテンシャル関数を用いなければ等電位線と電気力線の同時高速作画は出来ないのでは
ないかと思い、その関数を利用する方法ばかり考えていました。
今回ご教示頂いた事で、微係数を上手く利用すると複雑な作画も高速処理が出来るということを大変興味
深く勉強させていただきました。本当に有難う御座いました。

お礼の意味も込めて山中さんのコードを利用し纏めた全コードを記します。
!--1個の点電荷と接地導体球が存在する電場の電位・電気力線分布作図----
!<山中和義氏作の描画手法をまとめたもの>

!----描画座標・条件設定
LET Xa=-10
LET Xb=-Xa+4
LET Ya=Xa-2
LET Yb=-Ya+2
SET WINDOW Xa,Xb,Ya,Yb !描画範囲
DRAW grid
LET cEps=1E-3 !誤差精度設定

!----電荷位置設定と描画
LET r=1 !導体球の半径
LET a1=4 !点電荷のｘ座標
LET a2=r^2/a1 !影像電荷のx座標
LET q1=1 !点電荷の大きさ
LET q2=-r/a1*q1 !影像電荷の大きさ
DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数
SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)
SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

!----電気力線の描画----
SET LINE COLOR "black"
LET L=0.2 !移動量
LET h=0.001 !増分
FOR th=0 TO 360 STEP 3 !点電荷から放射状に
LET theta=RAD(th)
LET x=a1
LET y=0
PLOT LINES: x,y;
DO
WHEN EXCEPTION IN
LET w=f(x,y) !電位を得る
LET wx=f(x+h,y) !偏微分係数から
LET wy=f(x,y+h)
LET vx=-(wx-w)/h !∂f/∂x
LET vy=-(wy-w)/h !∂f/∂y
LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る
USE
END WHEN
LET x=x+L*COS(theta) !移動させる
LET y=y+L*SIN(theta)
PLOT LINES: x,y;
IF ABS(x)>Xb THEN EXIT DO !表示範囲外なら
IF ABS(y)>Yb THEN EXIT DO
IF ABS(x*x+y*y)<r THEN EXIT DO !接地導体球内なら
LOOP
PLOT LINES !次へ
NEXT th

!---画面仕上げ用（点電荷位置・接地導体球）再描画----
SET AREA COLOR "red" !点電荷の位置に赤丸印を表示
DRAW disk WITH SCALE(0.1)*SHIFT(a1,0)
SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示
DRAW disk WITH SCALE(r)*SHIFT(0,0)

!SET WINDOW -6,6,-6,6 !表示領域
!SET WINDOW -15,15,-15,15 !表示領域
!DRAW grid !座標

　　　　　　　　　　└プログラムの続き 島村1243 2007/12/16 19:05:45  ツリーへ

Re: 完璧に目的達成です！ 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/12/16 19:05:45

プログラムの続き

!---等電位分布作図(関数f(x,y)の等高線)-----
SET LINE COLOR "red"
LET h=0.001 !増分（偏微分係数）
LET d=0.015 !増分
FOR fc=0 TO 0.5 STEP d !等高線の関数値
CALL contour(fc,1,1,d)
NEXT fc

FOR fc=0 TO 0.06 STEP d !※左端
CALL contour(fc,-10,1,d)
NEXT fc

!--以下はSUB_routine------
SUB newton(fc,x,y, fx,fy,grad2) !ニュートン法
DO
LET ff=f(x,y) !∇f
LET fx=(f(x+h,y)-ff)/h
LET fy=(f(x,y+h)-ff)/h
LET grad2=fx*fx+fy*fy

IF grad2<1e-10 THEN
LET x=1e30
EXIT SUB
END IF

LET t=(fc-ff)/grad2
LET x=x+t*fx
LET y=y+t*fy
LOOP WHILE t*t*grad2>cEps*cEps
END SUB

SUB contour(fc,x,y,d) !等高線を描く
LET i=0
DO
CALL newton(fc,x,y, fx,fy,grad2)
IF ABS(x)+ABS(y)>1e10 THEN EXIT SUB
IF i=0 THEN
PLOT LINES: x,y; !始点
LET x0=x
LET y0=y
ELSE
PLOT LINES: x,y; !折れ線でつなげる
END IF
IF i>2 AND (x-x0)^2+(y-y0)^2<d*d THEN EXIT DO !始点近傍なら、終了

LET t=d/SQR(grad2)
LET x=x+fy*t
LET y=y-fx*t

LET i=i+1
LOOP
PLOT LINES: x0,y0 !閉じる
END SUB

END

インデックスへ

EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

もっと高速な方法はあるでしょうか？	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/13 14:21:55 ** この記事は1回修正されてます
電気影像法（映像法）の理屈と等角写像を組み合わせて、複数の点電荷が有る場合の２次元静電場の画像を描かせる方法を模索中です。 z平面の（a1,j0）と(a2,j0)の２点にそれぞれq1、q2の電荷（+,-を考慮）を置いた場合の等電位分布と電気力線分布図は、 w=q1Log(z-a1)+q2Log(z-a2)、w=u+jv （ただし左記以外の定数は省略）をz=f(w)と解析的に表す事が出来れば電位分布は j=sqr(-1) For u=-10 to 10 step 0.2 For v=-pi to pi step pi/20 w=u+jv z=f(w) plot lines:Re(z),Im(z); next plot lines next の様にして高速に、しかも綺麗に作図出来るはずです（電気力線分布も同様）。ちなみにa2=-a1,q1=q2=+1の場合（２個の正電荷がy軸対称に設置された場合）は w=Log(z-a1)+Log(z+a1) から z=f(w)=sqr(a1^2+exp(w))を得て作図が高速に行えました。しかしa1<>a2,q1<>q2の様な一般的な場合は解析解が得られない（と言うより解析方法が分からない）ので上記写像の理屈で作画が出来ません。電位分布だけなら各点電荷による電位を合成した式に、以前の投稿「根軌跡図の作成」で山中様に教えて頂いた作画方法 http://freebbs.around.ne.jp/article/b/basic/96/sfpfow/vzgpmu.html#vzgpmu を利用して下記の様に作成できます。

Re: もっと高速な方法はあるでしょうか？	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/13 14:23:05 ** この記事は1回修正されてます
続き !接地導体球と1個の点電荷による静電場の電位分布作図 !等式に基づくドット画像を描く山中氏の方法を利用 OPTION ARITHMETIC complex LET Xa=-6 LET Xb=-Xa LET Ya=-6 LET Yb=-Ya SET WINDOW Xa,Xb,Ya,Yb !描画範囲 DRAW grid ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算 LET c_Div=3 !分解数 LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る LET dy=(Yb-Ya)/b/c_Div LET cEps=1E-3 !誤差　※精度が気になる!? LET j=SQR(-1) !虚数単位 LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標セット LET a2=r^2/a1 !影像電荷のx座標計算 LET q1=1 !点電荷の大きさセット LET q2=-r/a1q1 !影像電荷の大きさ計算 !点電荷の位置に赤丸印を表示 SET AREA COLOR "red" DRAW disk WITH SCALE(0.1)SHIFT(a1,0) !接地導体球の位置に緑丸印を表示 SET AREA COLOR "green" DRAW disk WITH SCALE(r)SHIFT(0,0) SET POINT STYLE 1 FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定 IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red" END IF FOR k=0 TO c_Divb-1 !Ｙ座標　※ドットをもとに LET y=Ya+dyk ![Ya,Yb] FOR i=0 TO c_Diva-1 !Ｘ座標 LET x=Xa+dxi ![Xa,Xb] LET z=x+yj !空間上の点を計算する IF z<> a1 AND z<>a2 THEN LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位 IF ABS(w-wb) < cEps THEN PLOT POINTS: x,y !点画 END IF NEXT i NEXT k NEXT wb END しかし、（１）1.7GHzのPCでは作図に長い時間を要する。（２）ドット表示なので線が粗い。また、誤差1/1000設定でもドット幅が広い部分が生じてしまう。（３）電気力線分布が得られない。と言う難点があります。何か良い解決策は有るでしょうか？

Re: 続き	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/14 08:55:07
(1)について（処理手順の改良）等高線の数だけ複素平面を走査して毎回wを算出すると時間を費やします。 FOR文のネストをかえると多少速くなります。「１回の走査で該当する等高線を描く」ように改修する場合（FOR文の部分抜粋） SET POINT STYLE 1 FOR k=0 TO c_Divb-1 !Ｙ座標　※ドットをもとに LET y=Ya+dyk ![Ya,Yb] FOR i=0 TO c_Diva-1 !Ｘ座標 LET x=Xa+dxi ![Xa,Xb] LET z=x+y*j !空間上の点を計算する IF z<>a1 AND z<>a2 THEN LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位 FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定 IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red" END IF IF ABS(w-wb) < cEps THEN PLOT POINTS: x,y !点画 NEXT wb END IF NEXT i NEXT k

Re: (1)について（処理手順の改良）	返事を書くノートメニュー
SECOND <jjqdmekgpt> 2007/12/14 21:34:04
!もったいない、整理すれば、４倍速。 OPTION ARITHMETIC complex LET Xa=-6 LET Xb=-Xa LET Ya=-6 LET Yb=-Ya SET WINDOW Xa,Xb,Ya,Yb !描画範囲 DRAW grid ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算 LET c_Div=3 !分解数 LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る LET dy=(Yb-Ya)/b/c_Div LET cEps=1E-3 !誤差　※精度が気になる!? ! LET j=SQR(-1) !虚数単位　　…使わない LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標セット LET a2=r^2/a1 !影像電荷のx座標計算 LET q1=1 !点電荷の大きさセット LET q2=-r/a1q1 !影像電荷の大きさ計算 !点電荷の位置に赤丸印を表示 SET AREA COLOR "red" DRAW disk WITH SCALE(0.1)SHIFT(a1,0) !接地導体球の位置に緑丸印を表示 SET AREA COLOR "green" DRAW disk WITH SCALE(r)*SHIFT(0,0) SET POINT STYLE 1 FOR y=Ya TO Yb STEP dy ![Ya,Yb] FOR x=Xa TO Xb STEP dx ![Xa,Xb] LET z=complex(x,y) !空間上の点z IF z<>a1 AND z<>a2 THEN LET w=q1/ABS(z-a1)+q2/ABS(z-a2) !点zにおける合成電位 FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定 IF ABS(w-wb) < cEps THEN IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red" PLOT POINTS: x,y !点画 END IF NEXT wb END IF NEXT x NEXT y END

Re: !もったいない、整理すれば、４倍速。	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/14 22:33:52
確かに速いです。 SECONDさん、いつもご教示有難う御座います。確かに速く、ストレスが感じなくなりますね。 SECONDさんのテクニックで４倍速になる理由が何なのかを知るために、山中さんのコードとの違いを見たのですが（１）z=u+j*vを使用せず、z=complex(x,y)を使用している。（２）画面走査に「FOR k=　　」のループを使用せず、刻み幅設定は「STEP」で指定し、FORループを２段に削減している。でした。今回の場合、上記（１）と（２）のどちらが大きな効力を発揮しているのでしょうか？

Re: 確かに速いです。	返事を書くノートメニュー
SECOND <jjqdmekgpt> 2007/12/14 22:50:11 ** この記事は1回修正されてます
どっちでしょう？実行頻度の最も多い　最も内側のステートメントの処理を軽くする努力が、一番効果します。文の順番も重要です。それから、使用する変数の種類を減らす工夫をすると、不思議な事に、自然に文の方が誘導してくれる傾向があります。これは、山中氏のアイデアがよかったから、マグレですね。(^^;) それより、逆写像の式を求める方が、大事ですね。私にもわかりません。

Re: (1)について（処理手順の改良）	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/14 22:06:46
２０秒速くなりました。山中さん、ご教示有難う御座います。同じ全画面走査を何回も繰り返すよりも、全画面走査は1回の方が速くなると言う論理ですね。私のPCでは原作が1分５５秒かかって完了に対し、山中さんの改修コード適用の場合は1分３５秒でした。

Re: 続き	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/15 09:43:05
(3)について　各点での流線 !複素ポテンシャル OPTION ARITHMETIC complex LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(z)=q1LOG(z-a1)+q2LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み） SET WINDOW -6,6,-6,6 !描画範囲 DRAW grid SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)SHIFT(0,0) LET h=0.01 !増分 FOR x=-6 TO 6 STEP 0.5 !複素平面を走査する FOR y=-6 TO 6 STEP 0.5 WHEN EXCEPTION IN LET w=f(complex(x,y)) !ベクトル成分を得る　※流線 LET wx=f(complex(x+h,y)) !偏微分係数から LET wy=f(complex(x,y+h)) LET vy=-Im(wx-w)/h !∂ψ/∂x LET vx= Im(wy-w)/h !∂ψ/∂y LET ww=complex(vx,vy) !成分に応じたベクトルを得る DRAW arrow WITH ROTATE(arg(ww))SCALE(ABS(ww)0.7)SHIFT(x,y) !※0.7 調整要 USE END WHEN NEXT y NEXT x PICTURE arrow !矢印を描く PLOT LINES: 0,0; 1,0 PLOT LINES: 1,0; 0.8, 0.1 PLOT LINES: 1,0; 0.8,-0.1 END PICTURE END

Re: (3)について各点での流線	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/15 11:03:55
!等ポテンシャルを色の濃淡で表現する OPTION ARITHMETIC complex LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(z)=q1LOG(z-a1)+q2LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み） SET WINDOW -6,6,-6,6 !描画範囲 SET POINT STYLE 1 FOR x=-6 TO 6 STEP 0.02 !複素平面を走査する FOR y=-6 TO 6 STEP 0.02 WHEN EXCEPTION IN LET w=f(complex(x,y)) !LET c=MOD(Re(w)1000,256)/256 !電位　※1000 調整要 LET c=MOD(Im(w)1000,256)/256 !電気力線 SET COLOR MIX(1) 0,0,c !青の濃淡 PLOT POINTS: x,y USE END WHEN NEXT y NEXT x DRAW grid !座標 SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(0.1)*SHIFT(a2,0) END

Re: (3)について各点での流線	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/15 20:11:13
これを利用すると流れの方向がわかるので、タートル（亀）をそれに沿って移動させると、、、（川の流れに身をまかせる、、、感じです。） !タートルグラフィックスで電気力線を描く OPTION ARITHMETIC complex LET r=1 !導体球の半径 LET a1=4 !点電荷のx座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(z)=q1LOG(z-a1)+q2LOG(z-a2) !w=f(z)=φ+iψ　吹き出し（吸い込み） SET WINDOW -6,6,-6,6 !描画範囲 DRAW grid LET L=0.2 !移動量 LET h=0.001 !増分 FOR th=0 TO 360 STEP 5 !点電荷から放射状に LET theta=RAD(th) LET x=a1 LET y=0 PLOT LINES: x,y; DO WHEN EXCEPTION IN LET w=f(complex(x,y)) !ベクトル成分を得る　※流線 LET wx=f(complex(x+h,y)) !偏微分係数から LET wy=f(complex(x,y+h)) LET vy=-Im(wx-w)/h !∂ψ/∂x LET vx= Im(wy-w)/h !∂ψ/∂y LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る USE END WHEN LET x=x+LCOS(theta) !移動させる LET y=y+LSIN(theta) PLOT LINES: x,y; IF ABS(x)>6 THEN EXIT DO !表示範囲外なら IF ABS(y)>6 THEN EXIT DO IF ABS(xx+yy)<r THEN EXIT DO !接地導体球内なら LOOP PLOT LINES !次へ NEXT th SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)*SHIFT(0,0) END

Re: これを利用すると流れの方向がわかるので、...	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/15 23:15:50 ** この記事は1回修正されてます
とても数理的で良いですね。山中さん、いろいろなご検討有り難うございます。電位分布（または電気力線）を色の濃淡に変換する案は、何故、連続的に変化している電位に縞模様の濃淡が出きるのか？と言う点が理解できない状態ですが、視覚的にはとても素晴らしい発想ですね。プログラムコードを良く見て何故？を考えようと思います。最後に示して頂いたタートルグラフィックスは、数理的な扱いも分かり易く、かつ、当初求めていた描画の目的が達成出来る解決策だと思います。山中さんのプログラム例では、正電荷から負電荷に向かって電気力線が放出されると言う考え方で、１ヶ所からの電気力線放出のみを計算していました。正電荷が２個有る場合にも、同じ様な考え方で２ヶ所について放出計算をすれば正解が得られるかを試行してみます。さて、解決策を示して頂き胸わくわくだったのですが、パソコンの描画結果をよくよく眺めてみると、何かおかしいなと気づきました。１個の正点電荷とその近くに接地導体球が有る二次元空間の電位分布は、電気影像法に依れば正電荷とそれの影像電荷のそれぞれの電荷によって生じる電位分布を代数加算して求められます。私が最初に提起したプログラムはその方法を使って電位分布を計算していますが、実際に計算される接地導体球(半径１としました）の表面電位はゼロになり、かつ、接地導体球の周囲には球（完全球ではない）状の電位分布が描かれているので正しい結果を表している事が分かります。ところが、複素ポテンシャル関数w=q1Log(z-a1)+q2Log(z-a2)を使って山中さんが作成してくださった電位分布状況図について接地導体球（原点を中心とした半径１）の周辺について見ると、等電位線が接地導体球表面と交わる形に描かれています。これは大きさの有る接地導体球にかんしては正しい作図を表していない（点電荷の電位分布図としては正しい作図の様に思えますが）と言う事です。つまり、「接地導体球内に仮想された影像電荷と真の点電荷の全体の複素ポテンシャル関数をそれぞれの複素ポテンシャル関数を単純に加算（例えばq1Log(z-a1)+q2Log(z-a2)の様に）して求め、その式から電位分布・電気力線分布を作図すれば良い」と言う私の考え方は根本的に正しく無かったのではないか？と言うことです。

Re: とても数理的で良いですね。	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/16 07:37:08
接地導体球の作画に問題が有る様複素ポテンシャル関数はLog関数、一方原点を中心に半径１の接地導体円を描いた座標系はLogとは無関係。その辺に誤りの根本が有るように感じています。

Re: 接地導体球の作画に問題が有る様	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/16 10:44:04 ** この記事は1回修正されてます
対策案できました。山中さんがご教示くださったタートルグラフィックスの手法と、SECONDさんがご教示くださった高速描画手法を組み合わせて、（複素ポテンシャル関数は諦めて）普通のポテンシャル関数を使って描画したら、（電位分布は誤差の為に線に幅が生じていますが）正しい分布状態の作図に成功しましたので以下に示します。お二方には感謝申し上げます。 !----点電荷と接地導体球の電場（電位分布と電気力線分布）作図---- !-------山中和義氏とSECOND氏の方法を合体---- LET Xa=-10 LET Xb=-Xa+4 LET Ya=Xa-2 LET Yb=-Ya+2 SET WINDOW Xa,Xb,Ya,Yb !描画範囲 DRAW grid ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算 LET c_Div=3 !分解数 LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る LET dy=(Yb-Ya)/b/c_Div LET cEps=1E-3 !誤差　※精度が気になる!? LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数 SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)SHIFT(0,0) !--電位分布の作図--- SET POINT STYLE 1 FOR y=Ya TO Yb STEP dy ![Ya,Yb] FOR x=Xa TO Xb STEP dx ![Xa,Xb] IF x<>a1 AND x<>a2 THEN LET w=f(x,y) !点(x,y)における合成電位 FOR wb=0 TO 0.5 STEP 0.02 !等電位幅値の設定 IF ABS(w-wb) < cEps THEN IF wb<=0 THEN SET POINT COLOR "black" ELSE SET POINT COLOR "red" PLOT POINTS: x,y !点画 END IF NEXT wb END IF NEXT x NEXT y !----電気力線の描画---- LET L=0.2 !移動量 LET h=0.001 !増分 FOR th=0 TO 360 STEP 3 !点電荷から放射状に LET theta=RAD(th) LET x=a1 LET y=0 PLOT LINES: x,y; DO WHEN EXCEPTION IN LET w=f(x,y) !電位を得る LET wx=f(x+h,y) !偏微分係数から LET wy=f(x,y+h) LET vx=-(wx-w)/h !∂f/∂x LET vy=-(wy-w)/h !∂f/∂y LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る USE END WHEN LET x=x+LCOS(theta) !移動させる LET y=y+LSIN(theta) PLOT LINES: x,y; IF ABS(x)>Xb THEN EXIT DO !表示範囲外なら IF ABS(y)>Yb THEN EXIT DO IF ABS(xx+yy)<r THEN EXIT DO !接地導体球内なら LOOP PLOT LINES !次へ NEXT th !---画面仕上げ用（点電荷位置・接地導体球）再描画---- SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)*SHIFT(0,0) END

Re: 対策案できました。	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/16 13:56:21
ポテンシャル関数f(x,y)なら、電位分布は等高線アルゴリズムで描けます。 !関数f(x,y)の等高線を描く LET cEps=0.001 !許容誤差 LET h=0.001 !増分（偏微分係数） SUB newton(fc,x,y, fx,fy,grad2) !ニュートン法 DO LET ff=f(x,y) !∇f LET fx=(f(x+h,y)-ff)/h LET fy=(f(x,y+h)-ff)/h LET grad2=fxfx+fyfy IF grad2<1e-10 THEN LET x=1e30 EXIT SUB END IF LET t=(fc-ff)/grad2 LET x=x+tfx LET y=y+tfy LOOP WHILE ttgrad2>cEpscEps END SUB SUB contour(fc,x,y,d) !等高線を描く LET i=0 DO CALL newton(fc,x,y, fx,fy,grad2) IF ABS(x)+ABS(y)>1e10 THEN EXIT SUB IF i=0 THEN PLOT LINES: x,y; !始点 LET x0=x LET y0=y ELSE PLOT LINES: x,y; !折れ線でつなげる END IF IF i>2 AND (x-x0)^2+(y-y0)^2<dd THEN EXIT DO !始点近傍なら、終了 LET t=d/SQR(grad2) LET x=x+fyt LET y=y-fxt LET i=i+1 LOOP PLOT LINES: x0,y0 !閉じる END SUB !---------- ここまでがサブルーチン !SET WINDOW -6,6,-6,6 !表示領域 SET WINDOW -15,15,-15,15 !表示領域 DRAW grid !座標 LET r=1 !導体球の半径 LET a1=4 !点電荷のx座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数 SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)*SHIFT(0,0) LET d=0.015 !増分 FOR fc=0 TO 0.5 STEP d !等高線の関数値 CALL contour(fc,1,1,d) NEXT fc FOR fc=0 TO 0.06 STEP d !※左端 CALL contour(fc,-10,1,d) NEXT fc END

Re: ポテンシャル関数f(x,y)なら、電位分布は等...	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/16 19:04:03 ** この記事は1回修正されてます
完璧に目的達成です！山中さん、幾度ものご教示有難う御座いました。当初は複素ポテンシャル関数を用いなければ等電位線と電気力線の同時高速作画は出来ないのではないかと思い、その関数を利用する方法ばかり考えていました。今回ご教示頂いた事で、微係数を上手く利用すると複雑な作画も高速処理が出来るということを大変興味深く勉強させていただきました。本当に有難う御座いました。お礼の意味も込めて山中さんのコードを利用し纏めた全コードを記します。 !--1個の点電荷と接地導体球が存在する電場の電位・電気力線分布作図---- !<山中和義氏作の描画手法をまとめたもの> !----描画座標・条件設定 LET Xa=-10 LET Xb=-Xa+4 LET Ya=Xa-2 LET Yb=-Ya+2 SET WINDOW Xa,Xb,Ya,Yb !描画範囲 DRAW grid LET cEps=1E-3 !誤差精度設定 !----電荷位置設定と描画 LET r=1 !導体球の半径 LET a1=4 !点電荷のｘ座標 LET a2=r^2/a1 !影像電荷のx座標 LET q1=1 !点電荷の大きさ LET q2=-r/a1q1 !影像電荷の大きさ DEF f(x,y)=q1/SQR((x-a1)^2+y^2)+q2/SQR((x-a2)^2+y^2) !合成電位関数 SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)SHIFT(0,0) !----電気力線の描画---- SET LINE COLOR "black" LET L=0.2 !移動量 LET h=0.001 !増分 FOR th=0 TO 360 STEP 3 !点電荷から放射状に LET theta=RAD(th) LET x=a1 LET y=0 PLOT LINES: x,y; DO WHEN EXCEPTION IN LET w=f(x,y) !電位を得る LET wx=f(x+h,y) !偏微分係数から LET wy=f(x,y+h) LET vx=-(wx-w)/h !∂f/∂x LET vy=-(wy-w)/h !∂f/∂y LET theta=ANGLE(vx,vy) !成分に応じたベクトルを得る USE END WHEN LET x=x+LCOS(theta) !移動させる LET y=y+LSIN(theta) PLOT LINES: x,y; IF ABS(x)>Xb THEN EXIT DO !表示範囲外なら IF ABS(y)>Yb THEN EXIT DO IF ABS(xx+yy)<r THEN EXIT DO !接地導体球内なら LOOP PLOT LINES !次へ NEXT th !---画面仕上げ用（点電荷位置・接地導体球）再描画---- SET AREA COLOR "red" !点電荷の位置に赤丸印を表示 DRAW disk WITH SCALE(0.1)SHIFT(a1,0) SET AREA COLOR "green" !接地導体球の位置に緑丸印を表示 DRAW disk WITH SCALE(r)*SHIFT(0,0) !SET WINDOW -6,6,-6,6 !表示領域 !SET WINDOW -15,15,-15,15 !表示領域 !DRAW grid !座標

Re: 完璧に目的達成です！	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/12/16 19:05:45
プログラムの続き !---等電位分布作図(関数f(x,y)の等高線)----- SET LINE COLOR "red" LET h=0.001 !増分（偏微分係数） LET d=0.015 !増分 FOR fc=0 TO 0.5 STEP d !等高線の関数値 CALL contour(fc,1,1,d) NEXT fc FOR fc=0 TO 0.06 STEP d !※左端 CALL contour(fc,-10,1,d) NEXT fc !--以下はSUB_routine------ SUB newton(fc,x,y, fx,fy,grad2) !ニュートン法 DO LET ff=f(x,y) !∇f LET fx=(f(x+h,y)-ff)/h LET fy=(f(x,y+h)-ff)/h LET grad2=fxfx+fyfy IF grad2<1e-10 THEN LET x=1e30 EXIT SUB END IF LET t=(fc-ff)/grad2 LET x=x+tfx LET y=y+tfy LOOP WHILE ttgrad2>cEpscEps END SUB SUB contour(fc,x,y,d) !等高線を描く LET i=0 DO CALL newton(fc,x,y, fx,fy,grad2) IF ABS(x)+ABS(y)>1e10 THEN EXIT SUB IF i=0 THEN PLOT LINES: x,y; !始点 LET x0=x LET y0=y ELSE PLOT LINES: x,y; !折れ線でつなげる END IF IF i>2 AND (x-x0)^2+(y-y0)^2<dd THEN EXIT DO !始点近傍なら、終了 LET t=d/SQR(grad2) LET x=x+fyt LET y=y-fxt LET i=i+1 LOOP PLOT LINES: x0,y0 !閉じる END SUB END