with(plots): with(plottools): com:=(x,y)->x+I*y; m:=2; Z[1]:=com(1,1); Z[2]:=com(3,2); q[1]:=3; q[2]:=5; phi := (z) -> sum(q[i]*log(z-Z[i]), i=1..m); # contourplot(Re(phi(com(x,y))),x=0..4,y=0..4,grid=[50,50]); # contourplot(Re(phi(com(x,y))),x=-10..15,y=-10..15,grid=[70,70],thickness=2,coloring=[red,blue]); # contourplot(Re(phi(com(x,y))),x=-10..15,y=-10..15,thickness=2,coloring=[red,blue]); me_init := proc(pmax,a,Q) local k; for k to pmax do a[k] := -sum(q[i]*Z[i]^k/k, i=1..m); od; Q := sum(q[i], i=1..m); end; me := (p,z) -> Q*log(z) + sum(a[k]/z^k, k=1..p); me_init(10,a,Q); gr:=160; contourplot(Re(phi(com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`exact potential`); contourplot(Re(me(0,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=0`); contourplot(Re(me(1,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=1`); contourplot(Re(me(2,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=2`); contourplot(Re(me(3,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=3`); contourplot(Re(me(10,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=10`);