Phase portraits of 2-dimensional systems

asked 2012-10-14 15:02:40 +0100

dxvxd
23 ●1 ●3 ●6

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

add a comment

3

answered 2012-10-14 16:54:37 +0100

achrzesz
2768 ●5 ●25 ●60

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

Comments

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

calc314 ( 2012-10-14 21:36:35 +0100 )edit

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

dxvxd ( 2012-10-15 00:20:50 +0100 )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 ( 2012-10-15 02:00:11 +0100 )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 ( 2012-10-15 03:35:53 +0100 )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 ( 2012-10-15 09:24:59 +0100 )edit

see more comments

2

answered 2012-11-07 07:41:25 +0100

Andrei Halanay
21 ●3

updated 2012-11-07 09:02:40 +0100

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

add a comment

1

answered 2012-10-15 10:55:45 +0100

calc314
4191 ●22 ●48 ●113

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

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 ( 2012-10-15 13:00:18 +0100 )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 ( 2012-10-15 14:44:19 +0100 )edit

add a comment

1

answered 2015-03-20 08:11:54 +0100

rws
1149 ●2 ●20 ●44 http://google.com/+Ral...

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

add a comment

1

answered 2013-12-06 19:35:25 +0100

Andrew F
11 ●1 ●1

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

add a comment

Phase portraits of 2-dimensional systems

5 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Phase portraits of 2-dimensional systems edit

5 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Phase portraits of 2-dimensional systems