SageMath を用いて電磁気学などで使われる概念についてプロットしてみます。物理学的な厳密性は投げ出しているので差し引いて読んでください。
まず、電荷をA~Dと名付け、位置を xs に、電荷の値を qs に格納しておきます。
A = vector([0,0]) B = vector([5,0]) C = vector([0,5]) D = vector([5,5]) xs = [A,B,C,D] qs = [1,-2,-3,2]
電荷だけをプロットしてみます。
charges = [] for x,q in zip(xs,qs): if q>0: charges.append(text("+",x,bounding_box={'boxstyle':'circle', 'fc':'w'})) else: charges.append(text("-",x,bounding_box={'boxstyle':'circle', 'fc':'w'})) sum(charges)
電荷が４つあるときの電場の様子をベクトル場で可視化してみます。
def g_vector_field(xs,qs,x_min,x_max,y_min,y_max): x,y=var('x,y') v = vector([x,y]) E = 0 Q = 1 for p,q in zip(xs,qs): E += Q*q*(v-p)/((v-p).norm())^3 E = E/E.norm() return plot_vector_field(E(x,y),(x,x_min,x_max),(y,y_min,y_max))
上で定義した関数では、ベクトルを正規化して描くようにしました。
sage: vf = g_vector_field(xs,qs,-8,13,-8,13) sage: vf Launched png viewer for Graphics object consisting of 1 graphics primitive
次にさきほど定義した電荷のプロットと重ねてみます。
sage: vf+sum(charges) Launched png viewer for Graphics object consisting of 5 graphics primitives
電気力線を描くために必要な関数を以下で定義します。厳密にいうと本数を指定して描画するようにしているので電気力線ではありませんが。
import sys radius = 0.1 def close(p,xs,qs,x_min,x_max,y_min,y_max): for x,q in zip(xs,qs): if q<0 and (p-x).norm() < radius: return True if p[0] < x_min or x_max < p[0]: return True if p[1] < y_min or y_max < p[1]: return True return False def points(theta,S,xs,qs,x_min,x_max,y_min,y_max): Q=1 P = S + 0.1*vector([cos(theta),sin(theta)]) N = 1000 l = [S,P] for i in range(N): E = 0 for x,q in zip(xs,qs): E += Q*q*(P-x)/((P-x).norm())^3 E = E/E.norm() P += 0.1*vector([float(E[0]),float(E[1])]) l.append(P) sys.stdout.write("\r{}".format(i+1)) sys.stdout.flush() if close(P,xs,qs,x_min,x_max,y_min,y_max): break print("") return l def charge_plot(xs,qs): g = [] for x,q in zip(xs,qs): if q>0: g.append(text("+",x,bounding_box={'boxstyle':'circle', 'fc':'w'})) else: g.append(text("-",x,bounding_box={'boxstyle':'circle', 'fc':'w'})) return sum(g) def ef_line(num_line,xs,qs,cs,x_min,x_max,y_min,y_max): gs = [] for x,q in zip(xs,qs): if q > 0: gs.append([points(2*i*pi/num_line,x,xs,qs,x_min,x_max,y_min,y_max) for i in range(num_line)]) g = sum([sum([line(gs[i][j],xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,color=cs[i]) for j in range(num_line)]) for i in range(len(gs))]) return g + charge_plot(xs,qs)
正の電荷は A と D なので、この2点から出る電気力線を描いてみます。
sage: g_A = sum([line(points(i*pi/4,A,xs,qs,-8,13,-8,13)) for i in range(8)]) sage: g_D = sum([line(points(i*pi/4,D,xs,qs,-8,13,-8,13)) for i in range(8)]) sage: g_A+g_D+charge_plot(xs,qs) Launched png viewer for Graphics object consisting of 20 graphics primitives
さきほど可視化したベクトル場と重ねます。
sage: vf+g_A+g_D+charge_plot(xs,qs) Launched png viewer for Graphics object consisting of 21 graphics primitives
当たり前ですが、線と各点のベクトルが同じ方向を向いています。
最後に、D の電荷を +2 から +19 まで変化させたときの電気力線の変化を調べてみます。
sage: frames = [ef_line(24,xs,[1,-2,-3,i],["blue","red","green","purple"],-8,13,-8,13) for i in sxrange(2,20)]
+2 のとき
sage: frames[0] Launched png viewer for Graphics object consisting of 52 graphics primitives
+19 のとき
sage: frames[-1] Launched png viewer for Graphics object consisting of 52 graphics primitives
アニメーションにする。
sage: animate(frames) Launched gif viewer for Animation with 18 frames
