Creating a Makie Streamplot using the output from DifferentialEquations.jl

General Usage
, ,
1

So I was looking to plot the phase diagram for simple ode systems using Makie. Now in the past I have used the Streamplot to display simple phase diagrams where I had an analytic solution for the ODE, as in this example.

However, in most cases ODEs are solved numerically. So I wanted to see if I can use a Makie streamplot to plot the phase diagram for the numerical solutions to an ODE. I found this example of doing the process manually, using regular plots.

# Simple Pendulum Problem
using DifferentialEquations, Plots
# Constants
const g = 9.81
const L = 1.0

#Initial Conditions
u₀ = [0,π/2]
tspan = (0.0,6.3)

#Define the problem
function simplependulum(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

#Pass to solvers
prob = ODEProblem(simplependulum,u₀, tspan)
sol = solve(prob,Tsit5())

#Plot
plot(sol,linewidth=2,
     title ="Simple Pendulum Problem", 
     xaxis = "Time", 
     yaxis = "Height", 
     label = ["Theta","dTheta"])

p = plot(sol, vars = (1,2), xlims = (-9,9), 
         title = "Phase Space Plot", 
         xaxis = "Velocity", 
         yaxis = "Position", 
         leg=false)
function phase_plot(prob, u0, p, tspan=2pi)
    _prob = ODEProblem(prob.f,u0,(0.0,tspan))
    sol = solve(_prob, Vern9()) # Use Vern9 solver for higher accuracy
    plot!(p, sol, vars = (1,2), xlims = nothing, ylims = nothing)
end
for i in -4pi:pi/2:4π
    for j in -4pi:pi/2:4π
        phase_plot(prob, [j,i], p)
    end
end
plot(p,xlims = (-9,9))

Now in this example, the authors just discretize some range of initial conditions, solve for the trajectory of a timespan, and then plot the output.

In the case of the Makie Streamplot, I can provide a range of initial conditions, but I don’t have a way to feed the solver like a timespan or such. For the Streamplot, it seems like I need to provide a value as Point2f() or Point3f(), but of course the solved numerical solution of the ODE is an array. So would I just return values of the second timestep in the solution array as a Point2f? I say the second timestep, because the first timestep is undoubtedly the initial value.

If anyone know a good recipe or example for this, please share the knowledge.

2

Streamplot doesn’t take a solution t -> F(t0, x0, t) but the vector field f(x) = d/dt F(t0, x, t) as input and uses a simple Euler scheme to numerically solve. Doesn’t that already work for your task?

3

Ahh okay. I must have read the docs wrong. So I don’t have to worry about the integration. I was thinking that I had to choose the solver etc. Okay, well let me give that a try. Yeah, I probably just got lost in notation and was not thinking clearly. Thanks @mschauer .

5

@mschauer Okay, I got it to work. Sorry, I had a typo in my first attempt. But the code works
now.

Here is the working code.

	struct SimplePendulum{T}
    	g::T
    	L::T
	end

	params = SimplePendulum(g, L)

	function h(x, P::SimplePendulum)
		Point2f(
			x[2],
			-(P.g/P.L)*sin(x[1])
		)
	end
	
	h(x) = h(x, params)

	fig, ax, pl = streamplot(h, -4π:0.001:4π, -4π:0.001:4π, colormap = :magma)
	fig

Here is the final image. Thanks again for setting me straight.

pendulum
pendulum1200×900 299 KB

3 Likes
6

Awesome!