# Trying to plot freezes Julia?

Hello,

I am trying to create a plot using Plots.jl and have run some code but it has already been over 15 minutes with no response from Julia. The screen looks like: Is this normal? It seems unusual to me. I am using Julia version 1.0 and have updated all packages.

My code (which has been checked by others in Using DifferentialEquations.jl for a simple food-chain ) is:

``````using DifferentialEquations

f = @ode_def_bare FoodChain begin
du1 = u1*(1-u1)-(1/e)*y2*x2*u1^(1+q)/(u012^(1+q)+u1^(q+1))*u2
du2 = y2*x2*u1^(1+q)/(u012^(1+q)+u1^(1+q))*u2-(1/e)*y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x2*u2
du3 = y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x3*u3
end e y2 x2 q u012 y3 x3 u023
u0 = [0.01,0.01,0.01]
tspan = (0.0,10000.0)
p = (1,2,3,4,5,6,7,8)
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,saveat=0.05)

using Plots
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)
``````

I am trying to create a graph series like this, in Julia. This plot was created in R and took much less time.

It is very probable that I have written something wrong in the plotting instructions. I am new to this and want to learn by doing, but if Julia does not generate plots I cannot experiment I was able to run some examples from http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-1-:-Solving-Scalar-Equations-1 but now Julia is not responding. What am I doing wrong in running Plot.jl?

I think you’re missing a ‘)’:

sol=solve(prob,saveat=0.05

Good eye! It was a copy/paste error, I ran it with the ‘)’ in Julia.

(Took some time to compile the DifferentialEquations package)

The layout should be a tuple:

plot(sol,xlabel = “Time” ,ylabel = “Density”, layout = (3,1))

I can only currently plot with UnicodePlots, but with that it works fine and displays a plot (whether the plot itself is correct I do not know).

You could try the GR backend with rasterized output, `Plots.gr(fmt=:png)`. This allows me to plot hundreds of thousands of points without performance issues.

Sorry, I don’t know what this means. I know GR is a package, and this suggests it works with Plots.jl?

Thank you, it produces a graph now, but alas, something else is amiss… Try the `plotdensity` kwarg: http://docs.juliadiffeq.org/latest/basics/plot.html#Density-1

`plot(sol, plotdensity=1000)`

`plot(sol,xlabel = "Time" ,ylabel = "Density", plotdensity=1000, layout = (3,1))`

Does not seem to make the graph change. I think the plot density should be 200 000 (tstart = 0, tmax = 10 000, increases in increments of 0.05) but this does not change the look of the graphs either.

Plots.jl supports several backends that do the actual plotting. Instead of:

``````using Plots
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)
``````

You can use the `gr` backend by doing this:

``````using Plots
Plots.gr(fmt=:png)
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)
``````

In my experience `gr` handles large number of points very fast, but maybe this is not your issue.

1 Like

It doesn’t change the look of the graphs because the additional 199,000 points aren’t visible due to limited screen resolution. Your screen is only a few thousand pixels wide, and plot density beyond that is overkill.

The solver’s behaving correctly–is the equation itself free of typos? Do you have the original R code handy?

The original R code is

``````library("deSolve")
library("EMD")

model.rma <-function(t, x, parms)
{
b1<-x
b2<-x
b3<-x

db1 <- b1*(1-b1)-(1/e)*y2*x2*b1^(1+q)/(b012^(1+q)+b1^(q+1))*b2
db2 <- y2*x2*b1^(1+q)/(b012^(1+q)+b1^(1+q))*b2-(1/e)*y3*x3*b2^(1+q)/(b023^(1+q)+b2^(1+q))*b3-x2*b2
db3 <- y3*x3*b2^(1+q)/(b023^(1+q)+b2^(1+q))*b3-x3*b3

return(list(c(db1,db2,db3)))
}

tmax <- 10000
b1.0 <- 0.01    # initial values
b2.0 <- 0.01
b3.0 <- 0.01

e <- 1          # assimilation efficiency
y2 <- 2.009     # maximum ingestion rate (relative to metabolic rate) of species 2
x2 <- 0.4       # metabolic rate of species 2
y3 <- 5         # maximum ingestion rate (relative to metabolic rate) of species 3
x3 <- 0.08      # metabolic rate of species 3
q <-  0         # Hill-Exponent = 1 + q
b012 <- 0.16129 # half-saturation density
b023 <- 0.5     # half-saturation density

#parms <-c(e,y2,x2,q,b012,y3,x3,b023)
xstart <-c(b1=b1.0, b2=b2.0, b3=b3.0)

times <-seq(0, tmax, 0.05)

out1 <-as.data.frame(lsoda(xstart, times, model.rma, parms))

par(mfrow=c(3,1))
par(las=1)
par(pty="m")  # maximises plotting area (instead of squares)

# plot time series
plot(
c(out1\$time,out1\$time,out1\$time),
c(out1\$b1,out1\$b2,out1\$b3),
xlab="time",
ylab="top predator density",
type="n"
)
lines(out1\$time,out1\$b3,col="red")

plot(
c(out1\$time,out1\$time,out1\$time),
c(out1\$b1,out1\$b2,out1\$b3),
xlab="time",
ylab="consumer density",
type="n"
)
lines(out1\$time,out1\$b2,col="orange")

plot(
c(out1\$time,out1\$time,out1\$time),
c(out1\$b1,out1\$b2,out1\$b3),
xlab="time",
ylab="prey density",
type="n"
)
lines(out1\$time,out1\$b1,col="darkgreen")
``````

Ah I see! This will be useful for me in the future once I start dealing with more data. Thank you.

It’s definitely a resolution/selective plotting issue with the plot output - when I plot only the first solution (`plot(sol[1,:],xlabel = "Time" ,ylabel = "Density")`) I get the reverse problem with the line being up instead of down below, i.e. seemingly only plotting local maxima instead of local minima. This might be a restriction of me plotting in the console though. Can you try only plotting one solution?

It looks like You didn’t use the same ODE parameters, so no wonder the plots look different. No need to set the timestep so small, either: DiffEq.jl does a good job with adaptive timestepping.

``````p = (1, 2.009, 0.4, 0, 0.16129, 5, 0.08, 0.5)
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob)
``````

2 Likes

That ought to do it! I suspected the data too now. Here’s some stats from before:

``````julia> mean(sol[1,:])
0.999537010691792

julia> median(sol[1,:])
0.9999999999999999
``````

and after your code:

``````julia> mean(sol[1,:])
0.6955485905767033

julia> median(sol[1,:])
0.9031611072330993
``````

Oh my. Yes you are correct, somehow parameters were switched. Thank you and sorry for making a silly mistake.