2

When solving differential equations and plotting the results, how can I increase the number of data points that are plotted? I have

using DifferentialEquations
using Plots

function lorenz(du,u,p,t)
    du[1] = 10.0*(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [5.0;0.0;0.0]
tspan = (0.0,100000.0)
prob = ODEProblem(lorenz,u0,tspan)

sol = solve(prob)

plot(sol, vars = 1, xlims = (10,100000),xscale =:log)

Specifically, when using something like PyPlots I can use:

x = linspace(0,100000,num=10000)

where I set num = 10000 which increases the number of samples and allows for a higher resolution of data points for longer integration timespans.

The obvious answer is to use PyPlots, but I'm not sure if I can even use PyPlots with the DifferentialEquations package. It would be nice to know how this is done for Plots. (As depending on the function, some plots come out very jagged).

Dontdownvote
  • 25
  • 3
  • I edited the answer after @ChrisRackauckas comment about the sampling part. Please see the last revision. – hckr Oct 29 '18 at 13:02

1 Answers1

1

Plots.jl is actually a wrapper around plotting backends like PyPlot.jl and GR.jl (which are also wrappers).

You can make Plots.jl use PyPlot by calling pyplot(), provided that you have PyPlot.jl installed.

using Plots
pyplot() # switch to PyPlot
plot(sol.t[2:end], sol[1,2:end], xlim =(10,100000), xscale=:log10)

Note that I started from the second index, because at the first index t is 0, which creates problem with log-scale (even if xlim is set). This uses every data point in the solution. If you are on the terminal this will open up PyPlot GUI, so you can zoom as you wish. If you are on Juno you can use gui(plot(...)) to open the window.

You can sample the system at arbitrary time points. (with the DifferentialEquations' interpolant sol(t))

t = range(10, stop = 100000, num=10000) # to reach t=100000 you need to adjust `maxiters` accordingly in your call to `solve`
samples = sol(t) # sample using interpolator
plot(samples.t, samples.u[1,:], xscale=:log10)

You could also use the recipe without log-scale and with plotdensity option.

plot(sol, vars=1, plotdensity=1000)

See here for more examples: http://diffeq.sciml.ai/stable/basics/plot.html

Chris Rackauckas
  • 18,645
  • 3
  • 50
  • 81
hckr
  • 5,456
  • 1
  • 20
  • 31
  • 1
    For the PyPlot part, I would suggest using the interpolant `sol(t)`, where `t = range(0,stop=100000,length=10000)` as the real counterpart to what was asked for in the question. But yes, using the plot recipe with `plotdensity` is the preferred option. – Chris Rackauckas Oct 29 '18 at 02:12