Ask Your Question
2

Phase portraits of 2-dimensional systems

asked 2012-10-14 08:02:40 -0500

dxvxd gravatar image

I'm trying to plot solutions to two dimensional ordinary differential equations. I found that Sage makes it easy to plot a vector field and, using ode_solver(), I can plot solutions on top of the vector field by specifying an initial condition y_0 and some range of time to run (t_span).

However, this method I'm using seems to be quite ad hoc, as I have to choose the right initial conditions and time span / know a lot about my system in order to plot a nice picture. Let's make this more concrete:

Say I want to draw a nice phase portrait for

$\dot{x} = -y$

$\dot{y} = -x$

First I generate the vector field:

var('x y')
VF=plot_vector_field([-y,-x],[x,-2,2],[y,-2,2])

Then I use ode_solver() to plot solutions with initial conditions going around the unit circle:

T = ode_solver()
T.function=lambda t,y:[-y[1],-y[0]]
solutions = []
c = circle((0,0), 1, rgbcolor=(1,0,1))
for k in range(0,8):
    T.ode_solve(y_0=[cos(k*pi/4),sin(k*pi,t_span=[0,1],num_points=100)
    solutions.append(line([p[1] for p in T.solution]))

This generates the following picture:

But if I change run the system for one more unit of time (set t_span=[0,2]), the picture gets skewed:

This makes sense, of course, because the magnitude of the vectors along $y=-x$ get big as you get further away from the origin. Similarly, the trajectory along $y=x$ has trouble getting to the origin because the magnitude of those vectors get very small. Which all makes me think there's a better way to do this. Thoughts?

edit retag flag offensive close merge delete

5 answers

Sort by ยป oldest newest most voted
3

answered 2012-10-14 09:54:37 -0500

achrzesz gravatar image
sage: maxima('plotdf([-y,-x],[x,y],[x,-2,2],[y,-2,2])')

Every click in the obtained field gives you a new trajectory

Precise coordinates of an initial point can be provided in plot setup

New x y coordinates + Enter adds a new trajectory

edit flag offensive delete link more

Comments

This doesn't work for me in the notebook. How can I get it to work?

calc314 gravatar imagecalc314 ( 2012-10-14 14:36:35 -0500 )edit

@calc314: could you describe what happens when you enter the above code?

dxvxd gravatar imagedxvxd ( 2012-10-14 17:20:50 -0500 )edit

Hello, If I run the code with sage notebok running on my laptop, new window from wxmaxima program with phase portrait appears. Probably will not work over the web (like sagenb.org).

robert.marik gravatar imagerobert.marik ( 2012-10-14 19:00:11 -0500 )edit

sage: indicates that I'm using Sage terminal (in graphical environment, linux 64-bit machine, sage 5.3) After copy & paste (without "sage:" part) into notebook my answer also works for me A new window with a vector field appears. Every click in this new window gives me a new trajectory Using the second position of the tools in this window I can introduce new x y coordinates (two floats separated by a space) When I press Enter after introducing them I'm obtaining the corresponding trajectory It can be repeated many times If only one trajectory is needed then the initial point coordinates can be defined in the command: sage: maxima('plotdf([-y,-x],[x,y],[x,-2,2],[y,-2,2],[trajectory_at,1,0])')

achrzesz gravatar imageachrzesz ( 2012-10-14 20:35:53 -0500 )edit

I was getting a response of `1` from sagenb. But, I just tried it on a local installation on my Mac, and a Tcl window did launch and responds to my mouse. Very Nice!

calc314 gravatar imagecalc314 ( 2012-10-15 02:24:59 -0500 )edit
2

answered 2012-11-07 00:41:25 -0500

Andrei Halanay gravatar image

updated 2012-11-07 02:02:40 -0500

The following function allows you to move the initial point and determines the "corners" of the picture with respect to the maximum an minimum values of the solutions.

@interact 
    def syst(x0=slider(-2,2,0.1,0.5),y0=slider(-2,2,0.1,0.5)):
        x,y=var('x y')
        vect=[-y,-x]
        sol=desolve_odeint(vect,[x0,y0],srange(-4,4,0.1),[x,y])
        xmin=sol[0,0]
        xmax=sol[0,0]
        ymin=sol[0,1]
        ymax=sol[0,1]
        for i in range(0,len(sol)):
            if sol[i,0] > xmax:
                xmax=sol[i,0]
            if sol[i,0]< xmin:
                xmin=sol[i,0]
        for i in range(0,len(sol)):
           if sol[i,1] > ymax:
               ymax=sol[i,0]
           if sol[i,1]< ymin:
               ymin=sol[i,0]
        p1=plot_vector_field((vect[0],vect[1]),(x,xmin-2,xmax+2),(y,ymin-2,ymax+2),plot_points=60)
        p=line(zip(sol[:,0],sol[:,1]))
        (p+p1).show(aspect_ratio=1)

Of course the solution using maxima is the nicest.

edit flag offensive delete link more
1

answered 2015-03-20 02:11:54 -0500

For completeness a bunch of ready-made solutions for ODE problems can be found at: http://wiki.sagemath.org/interact/diffeq

edit flag offensive delete link more
1

answered 2012-10-15 03:55:45 -0500

calc314 gravatar image

Beyond the idea in the other answer, I think you'll just have to keep up with the max and min values from the solutions in your loop and then use these to give you the necessary dimensions for your vector field plot. This is what I've done in the past. I don't know of any automated way to choose nice initial conditions and time spans to get great pictures.

edit flag offensive delete link more

Comments

How can I restrict the max/min values in the solutions? Is there a uniform way to do this or restrict the axes' range (e.g., force the plot onto [-2,2]x[-2,2] )

dxvxd gravatar imagedxvxd ( 2012-10-15 06:00:18 -0500 )edit
1

You should be able to use the `show` command with `xmin,xmax,ymin,ymax`. For example, based on your code above, this works for me: `q=sum(solutions)+VF+c` `q.show(xmin=-2,xmax=2,ymin=-2,ymax=2)`

calc314 gravatar imagecalc314 ( 2012-10-15 07:44:19 -0500 )edit
1

answered 2013-12-06 12:35:25 -0500

Andrew F gravatar image

One thing you can do is just make sure your vector function always returns the same length, that is, to wrap it in a call to a function that unitizes it (that's a word, right?)

def unit(v):
    mag = sqrt(v[0]*v[0]+v[1]*v[1])
    return (v[0]/mag, v[1]/mag)
VF = plot_vector_field(unit([-y,-x]),[x,-2,2],[y,-2,2])

That's about what I did, anyway.

edit flag offensive delete link more

Your Answer

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

Add Answer

Question Tools

3 followers

Stats

Asked: 2012-10-14 08:02:40 -0500

Seen: 2,531 times

Last updated: Mar 20 '15