Phase portraits of 2-dimensional systems

2

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?

add a comment

3

answered 12 years ago

achrzesz
2803 ●5 ●25 ●61

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

link

Comments

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

calc314 ( 12 years ago )

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

dxvxd ( 12 years ago )

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 ( 12 years ago )

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 ( 12 years ago )

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 ( 12 years ago )

see more comments

2

answered 12 years ago

Andrei Halanay
21 ●3

updated 12 years ago

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.

link

add a comment

1

answered 12 years ago

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.

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 ( 12 years ago )

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 ( 12 years ago )

add a comment

1

answered 10 years ago

rws
1149 ●2 ●21 ●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

link

add a comment

1

answered 11 years ago

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.

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 savecancel

5 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Phase portraits of 2-dimensional systems