|
|
後記のプログラムは、x,y,z ３次元空間中に点電荷と接地球導体が有る場合のz=0平面上の等電位曲線を描くもので、以前、この掲示板を通じて山中さんにご教示頂いたニュートン法のコードを使用しており正しく動作します。
しかしコードの一部が理解出来ず、お教え頂ければ有り難いです。
このプログラム中の「SUB newton(fc,x,y, fx,fy,grad2) !ニュートン法」に記述されている
LET t=(fc-ff)/grad2
LET x=x+t*fx
LET y=y+t*fy
のgrad2は、ヤコビ行列の逆行列作成で生じる行列式と推測されます。
電位ff=f(x,y)は２変数関数なので、ニュートン法でヤコビ行列を作成する為には、x,yに関するもう１個別の関数g(x,y)=kが必要と思うのですが、その関数はどのように考えればよいのでしょうか？
!**************************************************************
! 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)
!---等電位分布作図(関数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
|
|