# How to plot field lines in sage

I want to plot the field lines of an electric field, take for example that of an electric dipole. How can I do this in sage?

Edit: Perhaps my question was formulated too much in physicist slang... Mathematically I just want to plot the integral curves of a vector field for example: let $\vec{r}(x,y) = (x,y)$, $\vec{E_1} = \frac{\vec{r}}{|r|^{3}}$ and $\vec{E_2}(x,y) = \vec{E_1}(x-1,y)$

Then I want to plot the field lines of $\vec{E} = \vec{E_1} - \vec{E_2}$

The field line plot should look similar to this one:

I think in mathematica this is called a stream plot: http://reference.wolfram.com/mathemat...

edit retag close merge delete

Sort by » oldest newest most voted

Thanks for the edit. See if something like this does what you want.

var('x y')
E1=vector([x,y])/(sqrt(x^2+y^2))^3
E2=E1.subs(x=x-1)
F=(E1-E2)

tstop=0.4
step=0.005
int_conds=[[0.5,0.1*k] for k in [1..3,-3..-1]]
solns=sum([line(desolve_odeint(F,init,srange(0,tstop,step),[x,y]),thickness=2) for init in int_conds])
solns2=sum([line(desolve_odeint(F,init,srange(0,-1*tstop,-1*step),[x,y]),thickness=2) for init in int_conds])

F=(E1-E2)/norm(E1-E2)
solns+solns2+plot_vector_field(F,(x,-0.2,1.2),(y,-1,1))

more

Thanks, that's almost what I want. The missing point is that I want not only one or three lines but something like I added in the original post above.

( 2012-09-29 03:28:25 -0500 )edit

Sorry I had two typos in my Edit, which I fixed just now, see above. You should not normalize F! Perhaps you could include those changes in your examples. Thanks.

( 2012-09-29 05:21:04 -0500 )edit

I've updated my code a good bit to include your edits and some solving backwards and forwards from the initial condition.

( 2012-09-30 04:10:01 -0500 )edit

I don't think there's a way to automate this totally and avoid putting in all of the initial conditions. Here is another version.

var('x y')
E1=vector([x,3*y])/(sqrt(x^2+y^2))^(3/2)
E2=E1.subs(x=x-1)
F=(E1-E2)/norm(E1-E2)
int_conds=[[0.05,0.001],[0.05,0.005],[0.05,-0.001],[0.05,-0.005]]
solns=[desolve_odeint(F,init,srange(0,1.5,0.05),[x,y]) for init in int_conds]
nhalf=len(solns[0])/3
a=[arrow(solns[i][nhalf],solns[i][nhalf+1]) for i in range(0,len(solns))]
b=map(lambda x: line(x,thickness=2),solns)
show(sum(a)+sum(b)+point((0,0),size=50)+point((1,0),size=50)+plot_vector_field(F,(x,-0.1,1.1),(y,-.5,.5)),axes=false)


more

People interested in this question will want to perhaps help out in making our ticket for streamline plots a reality.

more

In the attached file one can see the maxima version (code in my previous answer)

sdfplot.pdf

more

I'm not sure whether you already have equations to use or not. How you plot these all depends, of course, on how the equations are given. But, here is something. Using information in section D.73.3 of this webpage on quantum mechanics, we need to plot ${(x^2+z^2)}^{(3/2)}=c x^2$ for various values of $c$. The following code produces a plot.

var('x z')
sum([implicit_plot((x^2+z^2)^(3/2)==c*x^2,(x,-5,5),(z,-5,5)) for c in [1,2,3,4]])


more

That are the equipotential surfaces and not the field lines. You can assume that I know an analytic expression for the electric field strength at every point.

( 2012-09-28 21:40:28 -0500 )edit

Can you post the expression?

( 2012-09-29 01:37:47 -0500 )edit

See my edit

( 2012-09-29 02:29:04 -0500 )edit