Ask Your Question
0

How to plot field lines in sage

asked 2012-09-28 11:23:45 -0500

sagefan gravatar image

updated 2012-09-29 09:36:13 -0500

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:

field line plot example

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

edit retag flag offensive close merge delete

5 answers

Sort by ยป oldest newest most voted
1

answered 2012-09-29 02:59:39 -0500

calc314 gravatar image

updated 2012-09-30 04:09:15 -0500

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))
edit flag offensive delete link more

Comments

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.

sagefan gravatar imagesagefan ( 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.

sagefan gravatar imagesagefan ( 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.

calc314 gravatar imagecalc314 ( 2012-09-30 04:10:01 -0500 )edit
1

answered 2012-09-29 05:10:30 -0500

calc314 gravatar image

updated 2012-09-29 05:11:17 -0500

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)

image description

edit flag offensive delete link more
1

answered 2012-09-29 15:10:15 -0500

kcrisman gravatar image

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

edit flag offensive delete link more
1

answered 2012-09-29 12:24:16 -0500

achrzesz gravatar image

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

sdfplot.pdf

edit flag offensive delete link more
1

answered 2012-09-28 11:53:20 -0500

calc314 gravatar image

updated 2012-09-28 11:55:22 -0500

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]])

image description

edit flag offensive delete link more

Comments

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.

sagefan gravatar imagesagefan ( 2012-09-28 21:40:28 -0500 )edit

Can you post the expression?

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

See my edit

sagefan gravatar imagesagefan ( 2012-09-29 02:29:04 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-09-28 11:23:45 -0500

Seen: 1,571 times

Last updated: Sep 30 '12