Plotting a phase portrait of a differential equation

General Usage
,
1

Hey Folks.

I was just wondering if there was an existing package or tool for plotting the phase portrait of a system of differential equations, or a phase portrait plus direction field? I am thinking of a plot like this for the system of equations.

\frac{dx}{dt} = -x \\ \frac{dy}{dt} = 2y

The corresponding plot would be (taken from Haberman’s ODE textbook):

Selection_002

I did not see any mention of this in the plots documentation for the Julia DifferentialEquations package, but perhaps the plot is in a different package.

It is not too bad to develop a plot like this myself, but just was wondering if anyone already had developed this–since it is such a common plot for differential equations.

Note that there was a previous post on this topic, but I don’t believe it received a response.

5 Likes
2

You can use the vars keyword argument to plot in Plots, which will allow you to get a phase shift.

For the direction plots, you should probably roll one yourself…

3

It’s a common plot for 2D examples that are used in textbooks, but quite hard to visualize for a hundred dimensional ODE :slight_smile:. But I would love to have a recipe for this if someone makes a nice one. It probably needs Makie to do it well.

3 Likes
4

Haha, come on @ChrisRackauckas haven’t you seen those cool videos and visualizations that Robert Ghrist has been tweeting out lately. He has a bunch of phase portraits and orbit plots in those. Haha, you know he is a Julia user too–for his computational topology stuff he uses Eirene.jl :).

Makie is a good idea. I was thinking about visualizing some of these simple systems of equations, just to examine some economics and social science models. Those are simple enough to explore in just 3 dimensions. Plus, you gotta master the basics if you are going to do the cooler stuff.

1 Like
5

+1 here for a simple script/recipe for platting phase portraits. I’m interested to use this in a class.

3 Likes
6

I’d like to see something like this too, and if I get a chance to put something together useful using Make I’ll post it. Meanwhile I’m using the built-in MacOS Grapher app which seems to be little known but is very powerful (though frustratingly the save/load round trip seems to be very buggy – beware!). Nice video showing its use for phase portraits:

7

Nowadays, you can do something like:

using Makie 
using AbstractPlotting
using AbstractPlotting.MakieLayout
AbstractPlotting.inline!(true)
odeSol(x,y) = Point(-x, 2y) # x'(t) = -x, y'(t) = 2y
scene = Scene(resolution =(400,400))
streamplot!(scene, odeSol, -2..2, -2..2, colormap = :plasma, 
    gridsize= (32,32), arrow_size = 0.07)
save("odeField.png", scene)

odeField

14 Likes
8

And with little bit more of extra work you can even get something like this:

using Makie 
using AbstractPlotting
using AbstractPlotting.MakieLayout
AbstractPlotting.inline!(true)
testField(x,y) = Point(-x, 2y)
x =  -2:0.1:2
y =  -2:0.1:2
u(x,y) = -x
v(x,y) = 2y
z = [log10(sqrt(u(x,y)^2 + v(x,y)^2)) for x in x, y in y]

scene, layout = layoutscene(resolution = (500, 400))
ax = layout[1, 1] = LAxis(scene, aspect = DataAspect(), 
    xlabel = "u", ylabel = "v", xticklabelsize=14, yticklabelsize=14)
fs = heatmap!(ax, x, y, z, colormap = :magma)
streamplot!(ax, testField, x, y, colormap = :magma,
    gridsize= (32,32), arrow_size = 0.04)
limits!(ax, -2, 2, -2, 2)
cbar = layout[1,2] = LColorbar(scene, fs, label = "log10[(u²+v²)½]", 
    labelsize = 14, ticklabelsize=14)
cbar.width = 12
cbar.height = Relative(3.7/4)
scene
save("odeField2.png", scene)

odeField2

More examples are easy to generate. Take for instance the analysis of limit cycles’s stability

using Makie 
using AbstractPlotting
using AbstractPlotting.MakieLayout
AbstractPlotting.inline!(true)

nonStablePoincare(x,y) = Point(x*(x^2+y^2-1) - y*(x^2+y^2+1),y*(x^2+y^2-1) + x*(x^2+y^2+1))
stableVanDerPaul(x,y) = Point(y, (1-x^2)*y -x)
semiStable(x,y) = Point(-y+ x*(-1+x^2+y^2)^2, x+y*(-1+x^2+y^2)^2) 

scene, layout = layoutscene(resolution = (900, 400))
ax1 = layout[1, 1] = LAxis(scene, aspect = DataAspect(), 
    xlabel = "x", ylabel = "y", title = "non-stable")
ax2 = layout[1, 2] = LAxis(scene, aspect = DataAspect(), 
    xlabel = "x", title = "stable")
ax3 = layout[1, 3] = LAxis(scene, aspect = DataAspect(),
    xlabel = "x", title = "semi-stable")

streamplot!(ax1, nonStablePoincare, -4..4, -4..4, colormap = :plasma, 
    gridsize= (32,32), arrow_size = 0.07)

streamplot!(ax2, stableVanDerPaul, -4..4, -4..4, colormap = :viridis, 
    gridsize= (32,32), arrow_size = 0.07)
streamplot!(ax3, semiStable, -4..4, -4..4, colormap = :rainbow1, 
    gridsize= (32,32), arrow_size = 0.07)
hideydecorations!(ax2, grid = false)
hideydecorations!(ax3, grid = false)
limits!(ax1, -4, 4, -4, 4)
limits!(ax2, -4, 4, -4, 4)
limits!(ax3, -4, 4, -4, 4)
scene
save("odeField3.png", scene)

odeField3
odeField3900×400 148 KB

25 Likes
9

How about dy/dx =y+ x? I tried odeSol(x,y) = Point(1, tan(x + y)), but it doesn’t seem right.

Update: odeSol(x,y) = Point(1, x + y) worked perfectly, thanks!

10

Gorgeous! But your code is not working for me. First, it accuses that both AbstractPlots and Makie have the same function (such as Scene). I tried disambiguate it, without success.

1 Like
11

you need to go to the new version (new syntax):

https://lazarusa.github.io/BeautifulMakie/FlowFields/streamplots3/

5 Likes
12

Hey, can we make an animation such that if we put few dots on the field, we can see how their position change with time?

1 Like