Formatting Plots

Chris Rackauckas

Since the plotting functionality is implemented as a recipe to Plots.jl, all of the options open to Plots.jl can be used in our plots. In addition, there are special features specifically for differential equation plots. This tutorial will teach some of the most commonly used options. Let's first get the solution to some ODE. Here I will use one of the Lorenz ordinary differential equation. As with all commands in DifferentialEquations.jl, I got a plot of the solution by calling solve on the problem, and plot on the solution:

using DifferentialEquations, Plots, ParameterizedFunctions
gr()
lorenz = @ode_def Lorenz begin
  dx = σ*(y-x)
  dy = ρ*x-y-x*z
  dz = x*y-β*z
end σ β ρ

p = [10.0,8/3,28]
u0 = [1., 5., 10.]
tspan = (0., 100.)
prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 1345-element Array{Float64,1}:
   0.0                
   0.0354861341350177 
   0.06066394416099348
   0.10188862127421744
   0.14484947449428065
   0.1983564366367698 
   0.2504990626839506 
   0.3056767768177142 
   0.3545280034970155 
   0.40770977583939344
   ⋮                  
  99.43850921240367   
  99.50820376840878   
  99.5966810633544    
  99.68534828643361   
  99.77443728645414   
  99.84980869692284   
  99.9110153350651    
  99.96735878458976   
 100.0                
u: 1345-element Array{Array{Float64,1},1}:
 [1.0, 5.0, 10.0]              
 [2.31565, 5.89756, 9.40679]   
 [3.23779, 7.04103, 9.23368]   
 [4.99386, 9.83293, 9.62611]   
 [7.42116, 13.9492, 11.5823]   
 [11.4597, 19.7531, 18.1042]   
 [15.4761, 21.5109, 29.8871]   
 [16.4475, 13.1242, 40.9711]   
 [12.8778, 2.61892, 41.2525]   
 [7.13698, -3.09341, 35.5052]  
 ⋮                             
 [-0.565857, -0.921084, 12.316]
 [-0.938645, -1.68821, 10.2883]
 [-1.99467, -3.77821, 8.43264] 
 [-4.49924, -8.6807, 8.21365]  
 [-10.002, -18.2144, 14.267]   
 [-16.1431, -22.0169, 31.2682] 
 [-16.359, -10.5027, 42.8011]  
 [-10.8707, 0.963967, 39.9983] 
 [-7.09417, 3.84177, 35.957]
plot(sol)

Now let's change it to a phase plot. As discussed in the plot functions page, we can use the vars command to choose the variables to plot. Let's plot variable x vs variable y vs variable z:

plot(sol,vars=(1, 2, 3))

We can also choose to plot the timeseries for a single variable:

plot(sol,vars=[:x])

Notice that we were able to use the variable names because we had defined the problem with the macro. But in general, we can use the indices. The previous plots would be:

plot(sol,vars=(1,2,3))
plot(sol,vars=[1])

Common options are to add titles, axis, and labels. For example:

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
xaxis="Time (t)",yaxis="u(t) (in mm)",label=["X","Y","Z"])

Notice that series recipes apply to the solution type as well. For example, we can use a scatter plot on the timeseries:

scatter(sol,vars=[:x])

This shows that the recipe is using the interpolation to smooth the plot. It becomes abundantly clear when we turn it off using denseplot=false:

plot(sol,vars=(1,2,3),denseplot=false)

When this is done, only the values the timestep hits are plotted. Using the interpolation usually results in a much nicer looking plot so it's recommended, and since the interpolations have similar orders to the numerical methods, their results are trustworthy on the full interval. We can control the number of points used in the interpolation's plot using the plotdensity command:

plot(sol,vars=(1,2,3),plotdensity=100)

That's plotting the entire solution using 100 points spaced evenly in time.

plot(sol,vars=(1,2,3),plotdensity=10000)

That's more like it! By default it uses 100*length(sol), where the length is the number of internal steps it had to take. This heuristic usually does well, but unusually difficult equations it can be relaxed (since it will take small steps), and for equations with events / discontinuities raising the plot density can help resolve the discontinuity.

Lastly notice that we can compose plots. Let's show where the 100 points are using a scatter plot:

plot(sol,vars=(1,2,3))
scatter!(sol,vars=(1,2,3),plotdensity=100)

We can instead work with an explicit plot object. This form can be better for building a complex plot in a loop.

p = plot(sol,vars=(1,2,3))
scatter!(p,sol,vars=(1,2,3),plotdensity=100)
title!("I added a title")