An Intro to DifferentialEquations.jl

Chris Rackauckas

Basic Introduction Via Ordinary Differential Equations

This notebook will get you started with DifferentialEquations.jl by introducing you to the functionality for solving ordinary differential equations (ODEs). The corresponding documentation page is the ODE tutorial. While some of the syntax may be different for other types of equations, the same general principles hold in each case. Our goal is to give a gentle and thorough introduction that highlights these principles in a way that will help you generalize what you have learned.

Background

If you are new to the study of differential equations, it can be helpful to do a quick background read on the definition of ordinary differential equations. We define an ordinary differential equation as an equation which describes the way that a variable $u$ changes, that is

\[ u' = f(u,p,t) \]

where $p$ are the parameters of the model, $t$ is the time variable, and $f$ is the nonlinear model of how $u$ changes. The initial value problem also includes the information about the starting value:

\[ u(t_0) = u_0 \]

Together, if you know the starting value and you know how the value will change with time, then you know what the value will be at any time point in the future. This is the intuitive definition of a differential equation.

First Model: Exponential Growth

Our first model will be the canonical exponential growth model. This model says that the rate of change is proportional to the current value, and is this:

\[ u' = au \]

where we have a starting value $u(0)=u_0$. Let's say we put 1 dollar into Bitcoin which is increasing at a rate of $98\%$ per year. Then calling now $t=0$ and measuring time in years, our model is:

\[ u' = 0.98u \]

and $u(0) = 1.0$. We encode this into Julia by noticing that, in this setup, we match the general form when

f(u,p,t) = 0.98u
f (generic function with 1 method)

with $ u_0 = 1.0 $. If we want to solve this model on a time span from t=0.0 to t=1.0, then we define an ODEProblem by specifying this function f, this initial condition u0, and this time span as follows:

using DifferentialEquations
f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 1.0

To solve our ODEProblem we use the command solve.

sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 5-element Array{Float64,1}:
 0.0                
 0.10042494449239292
 0.35218603951893646
 0.6934436028208104 
 1.0                
u: 5-element Array{Float64,1}:
 1.0               
 1.1034222047865465
 1.4121908848175448
 1.9730384275623003
 2.664456142481452

and that's it: we have succesfully solved our first ODE!

Analyzing the Solution

Of course, the solution type is not interesting in and of itself. We want to understand the solution! The documentation page which explains in detail the functions for analyzing the solution is the Solution Handling page. Here we will describe some of the basics. You can plot the solution using the plot recipe provided by Plots.jl:

using Plots; gr()
plot(sol)

From the picture we see that the solution is an exponential curve, which matches our intuition. As a plot recipe, we can annotate the result using any of the Plots.jl attributes. For example:

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false

Using the mutating plot! command we can add other pieces to our plot. For this ODE we know that the true solution is $u(t) = u_0 exp(at)$, so let's add some of the true solution to our plot:

plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

In the previous command I demonstrated sol.t, which grabs the array of time points that the solution was saved at:

sol.t
5-element Array{Float64,1}:
 0.0                
 0.10042494449239292
 0.35218603951893646
 0.6934436028208104 
 1.0

We can get the array of solution values using sol.u:

sol.u
5-element Array{Float64,1}:
 1.0               
 1.1034222047865465
 1.4121908848175448
 1.9730384275623003
 2.664456142481452

sol.u[i] is the value of the solution at time sol.t[i]. We can compute arrays of functions of the solution values using standard comprehensions, like:

[t+u for (u,t) in tuples(sol)]
5-element Array{Float64,1}:
 1.0               
 1.2038471492789395
 1.7643769243364813
 2.6664820303831105
 3.664456142481452

However, one interesting feature is that, by default, the solution is a continuous function. If we check the print out again:

sol
retcode: Success
Interpolation: Automatic order switching interpolation
t: 5-element Array{Float64,1}:
 0.0                
 0.10042494449239292
 0.35218603951893646
 0.6934436028208104 
 1.0                
u: 5-element Array{Float64,1}:
 1.0               
 1.1034222047865465
 1.4121908848175448
 1.9730384275623003
 2.664456142481452

you see that it says that the solution has a order changing interpolation. The default algorithm automatically switches between methods in order to handle all types of problems. For non-stiff equations (like the one we are solving), it is a continuous function of 4th order accuracy. We can call the solution as a function of time sol(t). For example, to get the value at t=0.45, we can use the command:

sol(0.45)
1.5542610480553116

Controlling the Solver

DifferentialEquations.jl has a common set of solver controls among its algorithms which can be found at the Common Solver Options page. We will detail some of the most widely used options.

The most useful options are the tolerances abstol and reltol. These tell the internal adaptive time stepping engine how precise of a solution you want. Generally, reltol is the relative accuracy while abstol is the accuracy when u is near zero. These tolerances are local tolerances and thus are not global guarantees. However, a good rule of thumb is that the total solution accuracy is 1-2 digits less than the relative tolerances. Thus for the defaults abstol=1e-6 and reltol=1e-3, you can expect a global accuracy of about 1-2 digits. If we want to get around 6 digits of accuracy, we can use the commands:

sol = solve(prob,abstol=1e-8,reltol=1e-8)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 9-element Array{Float64,1}:
 0.0                
 0.04127492324135852
 0.14679917846877366
 0.28631546412766684
 0.4381941361169628 
 0.6118924302028597 
 0.7985659100883337 
 0.9993516479536952 
 1.0                
u: 9-element Array{Float64,1}:
 1.0               
 1.0412786454705882
 1.1547261252949712
 1.3239095703537043
 1.5363819257509728
 1.8214895157178692
 2.1871396448296223
 2.662763824115295 
 2.664456241933517

Now we can see no visible difference against the true solution:

plot(sol)
plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

Notice that by decreasing the tolerance, the number of steps the solver had to take was 9 instead of the previous 5. There is a trade off between accuracy and speed, and it is up to you to determine what is the right balance for your problem.

Another common option is to use saveat to make the solver save at specific time points. For example, if we want the solution at an even grid of t=0.1k for integers k, we would use the command:

sol = solve(prob,saveat=0.1)
retcode: Success
Interpolation: 1st order linear
t: 11-element Array{Float64,1}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0
u: 11-element Array{Float64,1}:
 1.0               
 1.1029627851292922
 1.2165269512238264
 1.341783821227542 
 1.4799379510586077
 1.6323162070541606
 1.8003833264983586
 1.9857565541588764
 2.1902158127997704
 2.4157257420844966
 2.664456142481452

Notice that when saveat is used the continuous output variables are no longer saved and thus sol(t), the interpolation, is only first order. We can save at an uneven grid of points by passing a collection of values to saveat. For example:

sol = solve(prob,saveat=[0.2,0.7,0.9])
retcode: Success
Interpolation: 1st order linear
t: 3-element Array{Float64,1}:
 0.2
 0.7
 0.9
u: 3-element Array{Float64,1}:
 1.2165269512238264
 1.9857565541588764
 2.4157257420844966

If we need to reduce the amount of saving, we can also turn off the continuous output directly via dense=false:

sol = solve(prob,dense=false)
retcode: Success
Interpolation: 1st order linear
t: 5-element Array{Float64,1}:
 0.0                
 0.10042494449239292
 0.35218603951893646
 0.6934436028208104 
 1.0                
u: 5-element Array{Float64,1}:
 1.0               
 1.1034222047865465
 1.4121908848175448
 1.9730384275623003
 2.664456142481452

and to turn off all intermediate saving we can use save_everystep=false:

sol = solve(prob,save_everystep=false)
retcode: Success
Interpolation: 1st order linear
t: 2-element Array{Float64,1}:
 0.0
 1.0
u: 2-element Array{Float64,1}:
 1.0              
 2.664456142481452

If we want to solve and only save the final value, we can even set save_start=false.

sol = solve(prob,save_everystep=false,save_start = false)
retcode: Success
Interpolation: 1st order linear
t: 1-element Array{Float64,1}:
 1.0
u: 1-element Array{Float64,1}:
 2.664456142481452

Note that similarly on the other side there is save_end=false.

More advanced saving behaviors, such as saving functionals of the solution, are handled via the SavingCallback in the Callback Library which will be addressed later in the tutorial.

Choosing Solver Algorithms

There is no best algorithm for numerically solving a differential equation. When you call solve(prob), DifferentialEquations.jl makes a guess at a good algorithm for your problem, given the properties that you ask for (the tolerances, the saving information, etc.). However, in many cases you may want more direct control. A later notebook will help introduce the various algorithms in DifferentialEquations.jl, but for now let's introduce the syntax.

The most crucial determining factor in choosing a numerical method is the stiffness of the model. Stiffness is roughly characterized by a Jacobian f with large eigenvalues. That's quite mathematical, and we can think of it more intuitively: if you have big numbers in f (like parameters of order 1e5), then it's probably stiff. Or, as the creator of the MATLAB ODE Suite, Lawrence Shampine, likes to define it, if the standard algorithms are slow, then it's stiff. We will go into more depth about diagnosing stiffness in a later tutorial, but for now note that if you believe your model may be stiff, you can hint this to the algorithm chooser via alg_hints = [:stiff].

sol = solve(prob,alg_hints=[:stiff])
retcode: Success
Interpolation: specialized 3rd order "free" stiffness-aware interpolation
t: 8-element Array{Float64,1}:
 0.0                
 0.05653299582822294
 0.17270731152826024
 0.3164602871490142 
 0.5057500163821153 
 0.7292241858994543 
 0.9912975001018789 
 1.0                
u: 8-element Array{Float64,1}:
 1.0               
 1.0569657840332976
 1.1844199383303913
 1.3636037723365293
 1.6415399686182572
 2.043449143475479 
 2.6418256160577602
 2.6644526430553808

Stiff algorithms have to solve implicit equations and linear systems at each step so they should only be used when required.

If we want to choose an algorithm directly, you can pass the algorithm type after the problem as solve(prob,alg). For example, let's solve this problem using the Tsit5() algorithm, and just for show let's change the relative tolerance to 1e-6 at the same time:

sol = solve(prob,Tsit5(),reltol=1e-6)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 10-element Array{Float64,1}:
 0.0                 
 0.028970819746309166
 0.10049147151547619 
 0.19458908698515082 
 0.3071725081673423  
 0.43945421453622546 
 0.5883434923759523  
 0.7524873357619015  
 0.9293021330536031  
 1.0                 
u: 10-element Array{Float64,1}:
 1.0               
 1.0287982807225062
 1.1034941463604806
 1.2100931078233779
 1.351248605624241 
 1.538280340326815 
 1.7799346012651116
 2.0905717422346277
 2.4861021714470244
 2.6644562434913373

Systems of ODEs: The Lorenz Equation

Now let's move to a system of ODEs. The Lorenz equation is the famous "butterfly attractor" that spawned chaos theory. It is defined by the system of ODEs:

\[ \begin{align} \frac{dx}{dt} &= \sigma (y - x)\\ \frac{dy}{dt} &= x (\rho - z) -y\\ \frac{dz}{dt} &= xy - \beta z \end{align} \]

To define a system of differential equations in DifferentialEquations.jl, we define our f as a vector function with a vector initial condition. Thus, for the vector u = [x,y,z]', we have the derivative function:

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
lorenz! (generic function with 1 method)

Notice here we used the in-place format which writes the output to the preallocated vector du. For systems of equations the in-place format is faster. We use the initial condition $u_0 = [1.0,0.0,0.0]$ as follows:

u0 = [1.0,0.0,0.0]
3-element Array{Float64,1}:
 1.0
 0.0
 0.0

Lastly, for this model we made use of the parameters p. We need to set this value in the ODEProblem as well. For our model we want to solve using the parameters $\sigma = 10$, $\rho = 28$, and $\beta = 8/3$, and thus we build the parameter collection:

p = (10,28,8/3) # we could also make this an array, or any other type!
(10, 28, 2.6666666666666665)

Now we generate the ODEProblem type. In this case, since we have parameters, we add the parameter values to the end of the constructor call. Let's solve this on a time span of t=0 to t=100:

tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: [1.0, 0.0, 0.0]

Now, just as before, we solve the problem:

sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 1294-element Array{Float64,1}:
   0.0                  
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624077544510573
   0.009058075635317072 
   0.01695646895607931  
   0.0276899566248403   
   0.041856345938267966 
   0.06024040228733675  
   0.08368539694547242  
   ⋮                    
  99.39403070915297     
  99.47001147494375     
  99.54379656909015     
  99.614651558349       
  99.69093823148101     
  99.78733023233721     
  99.86114450046736     
  99.96115759510786     
 100.0                  
u: 1294-element Array{Array{Float64,1},1}:
 [1.0, 0.0, 0.0]                                                  
 [0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8]
 [0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6] 
 [0.9693591634199452, 0.08977060667778931, 0.0001438018342266937] 
 [0.9242043615038835, 0.24228912482984957, 0.0010461623302512404] 
 [0.8800455868998046, 0.43873645009348244, 0.0034242593451028745] 
 [0.8483309877783048, 0.69156288756671, 0.008487623500490047]     
 [0.8495036595681027, 1.0145425335433382, 0.01821208597613427]    
 [0.9139069079152129, 1.4425597546855036, 0.03669381053327124]    
 [1.0888636764765296, 2.052326153029042, 0.07402570506414284]     
 ⋮                                                                
 [12.999157033749652, 14.10699925404482, 31.74244844521858]       
 [11.646131422021162, 7.2855792145502845, 35.365000488215486]     
 [7.777555445486692, 2.5166095828739574, 32.030953593541675]      
 [4.739741627223412, 1.5919220588229062, 27.249779003951755]      
 [3.2351668945618774, 2.3121727966182695, 22.724936101772805]     
 [3.310411964698304, 4.28106626744641, 18.435441144016366]        
 [4.527117863517627, 6.895878639772805, 16.58544600757436]        
 [8.043672261487556, 12.711555298531689, 18.12537420595938]       
 [9.97537965430362, 15.143884806010783, 21.00643286956427]

The same solution handling features apply to this case. Thus sol.t stores the time points and sol.u is an array storing the solution at the corresponding time points.

However, there are a few extra features which are good to know when dealing with systems of equations. First of all, sol also acts like an array. sol[i] returns the solution at the ith time point.

sol.t[10],sol[10]
(0.08368539694547242, [1.0888636764765296, 2.052326153029042, 0.07402570506
414284])

Additionally, the solution acts like a matrix where sol[j,i] is the value of the jth variable at time i:

sol[2,10]
2.052326153029042

We can get a real matrix by performing a conversion:

A = Array(sol)
3×1294 Array{Float64,2}:
 1.0  0.999643     0.996105    0.969359     …   4.52712   8.04367   9.97538
 0.0  0.000998805  0.0109654   0.0897706        6.89588  12.7116   15.1439 
 0.0  1.78143e-8   2.14696e-6  0.000143802     16.5854   18.1254   21.0064

This is the same as sol, i.e. sol[i,j] = A[i,j], but now it's a true matrix. Plotting will by default show the time series for each variable:

plot(sol)

If we instead want to plot values against each other, we can use the vars command. Let's plot variable 1 against variable 2 against variable 3:

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

This is the classic Lorenz attractor plot, where the x axis is u[1], the y axis is u[2], and the z axis is u[3]. Note that the plot recipe by default uses the interpolation, but we can turn this off:

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

Yikes! This shows how calculating the continuous solution has saved a lot of computational effort by computing only a sparse solution and filling in the values! Note that in vars, 0=time, and thus we can plot the time series of a single component like:

plot(sol,vars=(0,2))

Internal Types

The last basic user-interface feature to explore is the choice of types. DifferentialEquations.jl respects your input types to determine the internal types that are used. Thus since in the previous cases, when we used Float64 values for the initial condition, this meant that the internal values would be solved using Float64. We made sure that time was specified via Float64 values, meaning that time steps would utilize 64-bit floats as well. But, by simply changing these types we can change what is used internally.

As a quick example, let's say we want to solve an ODE defined by a matrix. To do this, we can simply use a matrix as input.

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 11-element Array{Float64,1}:
 0.0                
 0.04315877410281615
 0.10749821507777302
 0.1798183032337637 
 0.27732622431309095
 0.3943300853390046 
 0.5271275340166153 
 0.6593505063074551 
 0.7998387276496801 
 0.9661257382228796 
 1.0                
u: 11-element Array{Array{Float64,2},1}:
 [0.306624934095197 0.5338066435533153; 0.19874765604983802 0.6868676641729
416; 0.6862329122068349 0.13858963914056321; 0.2031754356279607 0.805431744
4734056]   
 [0.2607854140452007 0.3613260168142186; 0.3061733660115813 0.6130920806716
098; 0.6484554422419746 0.09909176566962785; 0.33528190215003106 0.97147927
48298231]  
 [0.13596497914872052 0.02539146709353668; 0.39275624859948677 0.4140338403
5184776; 0.6235490412546256 0.1171227635342626; 0.5191531472173252 1.191169
1282105947]
 [-0.08267310627767702 -0.4573896273601639; 0.3952152446091268 0.0883641473
6267606; 0.6581128735546947 0.27062599647166163; 0.6975117123588294 1.38567
9442209257]
 [-0.49551765265574466 -1.2570950504304415; 0.27500883360229983 -0.45014272
42424515; 0.8439611836541472 0.7435773073443797; 0.8693546862420823 1.53201
7905622345]
 [-1.1251069448214652 -2.3625200085353635; 0.037325541951849256 -1.08123128
50121293; 1.3260226380299647 1.765626779016673; 0.9330006745006937 1.480132
4753451832]
 [-1.9045113695239797 -3.6235436930157006; -0.16220499769796015 -1.46363209
84855904; 2.246749925213862 3.5332898325322093; 0.7601909703615497 1.043513
9453705078]
 [-2.5636091056660977 -4.5900120441897645; -0.02547223902998061 -1.09470953
18805466; 3.510744970192141 5.809641306139268; 0.2774234486076798 0.1395762
174433236] 
 [-2.8588358769826323 -4.8675661833207755; 0.7924931322952213 0.55502678447
71868; 5.042857181653392 8.434977730075753; -0.6037853546431378 -1.36780693
54307625]  
 [-2.2056076203133 -3.4860281122307057; 3.002992564478892 4.630325770000123
; 6.583433020832067 10.916694329795456; -2.0967627480607884 -3.801465962226
624]       
 [-1.8905678647543251 -2.903621089761297; 3.6371785566754515 5.759661435694
562; 6.784809971693031 11.212684216367514; -2.446836669801699 -4.3599865292
66994]

There is no real difference from what we did before, but now in this case u0 is a 4x2 matrix. Because of that, the solution at each time point is matrix:

sol[3]
4×2 Array{Float64,2}:
 0.135965  0.0253915
 0.392756  0.414034 
 0.623549  0.117123 
 0.519153  1.19117

In DifferentialEquations.jl, you can use any type that defines +, -, *, /, and has an appropriate norm. For example, if we want arbitrary precision floating point numbers, we can change the input to be a matrix of BigFloat:

big_u0 = big.(u0)
4×2 Array{BigFloat,2}:
 0.306625  0.533807
 0.198748  0.686868
 0.686233  0.13859 
 0.203175  0.805432

and we can solve the ODEProblem with arbitrary precision numbers by using that initial condition:

prob = ODEProblem(f,big_u0,tspan)
sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 6-element Array{Float64,1}:
 0.0                
 0.11309299507432245
 0.34314692151136   
 0.638787912178119  
 0.9249800093396419 
 1.0                
u: 6-element Array{Array{BigFloat,2},1}:
 [0.306624934095196977068553678691387176513671875 0.53380664355331530757098
3350160531699657440185546875; 0.1987476560498380173669374926248565316200256
34765625 0.6868676641729416498804994262172840535640716552734375; 0.68623291
22068349107621543225832283496856689453125 0.1385896391405632144255832827184
3492984771728515625; 0.203175435627960698781180326477624475955963134765625 
0.805431744473405597517512433114461600780487060546875]                     
                                                                           
                                                                           
                                                               
 [0.12195360456574192164308358623413396755370929614133176343692674857595920
89283282 -0.008123130290881944410204034456520484439167574288627050549707058
882818710342921869; 0.39634232463333060475342878723397023611671692453560022
20324236305738155185910297 0.3922711469618407524571385921694548208680191059
529293489153371312310947138743918; 0.62360874142229614660269413816440946870
87370580824993887652107972174806719590678 0.1236409384743709259797972837855
357447990439632298887694654224199400105772431161; 0.53416347348269343867994
57018811739572433059005299242711070019148084150545073337 1.2083846743849025
18564024515351644107890176707970397695234248206414887624922368]
 [-0.8361838634484524348538168488462141801511078402805525410153586829197726
692808571 -1.86702009280275907517976463220718291159544885656355164272253725
3422383845174623; 0.1453811433020334574482407049325057368670470738112201107
235484165871856652496102 -0.82134478429799790787504938541531905469749406674
00386286734161179476112442392394; 1.077771068365623680961642145379083157281
57412606957580092742711657289164482844 1.2553877009621635276525343870881896
86801854167614311182221002263583654514473809; 0.927145045985413897479125565
3032764995985712313551139615976505988068257946842391 1.53723087864554990641
4797277498265984897737315208334480701715945474510318806641]    
 [-2.4781841695365935139189432377168542891767829957296609122941866219442182
78896451 -4.473442094117027675809815446775958785211969683715454078985814799
39549300369079; -0.08006364569705258042016347463273375867687218500322261325
409698105571951190333812 -1.21777113075278823531712886538386589247841244623
8940453398155952886774203517709; 3.2966315204870468195194565552806077969134
38470665102008643034908539507467653621 5.4321316110811692416608908089319752
02884242133404203141998361295001144659681931; 0.374349345424851256103484671
6131550426370925361654115708204641548117654790737735 0.31282775172215700341
93926449343762287138753819339916998491633864398598808831969]   
 [-2.4975818474654782027132192616937496987044671987510889796984588261379768
8739931 -4.0447100330112030751965498675045126574427038896461677180992339776
29276787366896; 2.317330136264621470670536733728552342915045303595855998302
506298055630934243305 3.393959417472666081039083941359688764565665061035505
721605949015682598273199206; 6.27290836063364715834215598699423777433297061
3699990873953739050279600565477863 10.4376619318291287301568866309489208220
0057863455139268312561945612660615303942; -1.689055448264079339133670789291
3297438964986327492426444109202910454414785522 -3.1459536741694805543316146
72864919967432018745488031128508039324887463085414649]         
 [-1.8905635653064640482754048965318043472960738675217969596201048475534655
62833882 -2.903613161693912596462335234152228047992884765645552376599867304
615746874344366; 3.63719468066918806341320932471563124940388939900708102008
3798980087192646198465 5.75968789300718216160605225298983954645231221977570
3753148103416886321220703276; 6.7848191643649295916304519941107790356057684
47698696199506860813363339546739066 11.212699088749110671833971366870494546
78391245938108692830749133496976192597072; -2.44684488609446281136992003045
2810812989399961122674021029527293419007432081828 -4.3600004123377005404380
03490097475473389116618818531156433722540001840908844477]
sol[1,3]
-0.836183863448452434853816848846214180151107840280552541015358682919772669
2808571

To really make use of this, we would want to change abstol and reltol to be small! Notice that the type for "time" is different than the type for the dependent variables, and this can be used to optimize the algorithm via keeping multiple precisions. We can convert time to be arbitrary precision as well by defining our time span with BigFloat variables:

prob = ODEProblem(f,big_u0,big.(tspan))
sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 6-element Array{BigFloat,1}:
 0.0                                                                       
       
 0.113092995074322453107164827202916505594705516625996223111945149913284021
3685782
 0.343146921511359988812756981258199791752946439061231935036208672177190552
5139586
 0.638787912178119045475686716403453696145703840976671770300362717636169010
3911218
 0.924980009339641787803025370802957820064305280645527971524647529343989238
2198667
 1.0                                                                       
       
u: 6-element Array{Array{BigFloat,2},1}:
 [0.306624934095196977068553678691387176513671875 0.53380664355331530757098
3350160531699657440185546875; 0.1987476560498380173669374926248565316200256
34765625 0.6868676641729416498804994262172840535640716552734375; 0.68623291
22068349107621543225832283496856689453125 0.1385896391405632144255832827184
3492984771728515625; 0.203175435627960698781180326477624475955963134765625 
0.805431744473405597517512433114461600780487060546875]                     
                                                                           
                                                                           
                                                             
 [0.12195360456574192057103450125829479835154212449537040631419147193227177
07072281 -0.008123130290881946954846439953346604247685795531381076954615265
216926057964851296; 0.39634232463333060500035154968630411236237809056373749
67291028820745967835867344 0.3922711469618407507967721492263185395342890073
81766330685625774594538169690652; 0.623608741422296146622188495992379377326
1012111860142087916384383100848754432546 0.12364093847437092650170876507509
88252011078298374602277522448929137058846502853; 0.534163473482693439801594
9598190797146120758990521065990311201727125884303541186 1.20838467438490251
984570602012818073090786122054696368793424589674326470922365]
 [-0.8361838634484524302515241377868298585150783921197594796200281898496975
666674019 -1.86702009280275906714481657557557772254317523145257578870427033
7325883704383484; 0.1453811433020334592194198867700128021358065155340114402
40125065565642004128116 -0.821344784297997903320173617608476524831684965801
178391546925474249475267299101; 1.07777106836562367736866494979442563453648
6095469728594798970600327565775724293 1.25538770096216352007838104928602742
9091119897403465064103255124923206814475632; 0.9271450459854138970877764411
093660000440062290278553621204170193599798676166655 1.537230878645549906894
162691393353603227149486073413270940164965502830507778154]   
 [-2.4781841695365933879454296869576851158794725492090979510087154475624816
40308525 -4.473442094117027500962149810242483127063290604571647473071064187
993636195979255; -0.0800636456970526473426825590296719382059392407749978855
3243695332637183329478918 -1.2177711307527883897259917368996975453640322015
70350741827901448233547018592178; 3.296631520487046521606572939309413704811
029124693202574218616309645150411280569 5.432131611081168714399861849631519
193841227727803785362948080853293664273521618; 0.37434934542485138684240887
51458359168829629652339837772895985107618118489826407 0.3128277517221572388
320863277424571593898950652132176086684904308392850230973752]
 [-2.4975818474654785620583333419028132327401174606322107203958802825370885
2923473 -4.0447100330112037811789973414846617842896582299344829841646119747
41954017475663; 2.317330136264620532155919713301594091246140146802376310965
55742941879474848447 3.3939594174726643759505719833718108393819666973928514
25070913540376507610733073; 6.272908360633646656799832472605245173020138797
91124780627884217093459523595315 10.437661931829127942742116637101708750848
09021359675225845126531385156749683626; -1.68905544826407875647187275626951
5308655502536388607395692873313231436862205155 -3.1459536741694796133885914
09304003319814855469322625609339840227466509054260104]       
 [-1.8905635653064639047282991343505179655840798093162889698003912501902176
05059431 -2.903613161693912334222297344693476087420852212725142231232173351
309153331868035; 3.63719468066918833601693286989159200951330508996604387659
1769828411954450994223 5.75968789300718264451102885315708689832163972296572
6432333575598555749412180328; 6.7848191643649296626210778813473088663956547
37119855069667803522335464117944441 11.212699088749110772509816303551421403
75641742621155666612423517643312409626039; -2.44684488609446295706038840307
4582141844252660843255139959116702928450319605342 -4.3600004123377007720866
8282413065848764416194625002658163522925887031917076989]

Let's end by showing a more complicated use of types. For small arrays, it's usually faster to do operations on static arrays via the package StaticArrays.jl. The syntax is similar to that of normal arrays, but for these special arrays we utilize the @SMatrix macro to indicate we want to create a static array.

using StaticArrays
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 11-element Array{Float64,1}:
 0.0                
 0.04373840556447778
 0.10504494381022488
 0.1832067917446593 
 0.27624508173685597
 0.38220516978577257
 0.5082351570065224 
 0.6464984228378687 
 0.7904627863642966 
 0.9604154264967177 
 1.0                
u: 11-element Array{StaticArrays.SArray{Tuple{4,2},Float64,2,8},1}:
 [0.9556033595934601 0.4445826506822814; 0.6653733603188241 0.3599066578409
216; 0.08928453186876584 0.9161879858955384; 0.5023138864356134 0.410089319
6093318]    
 [0.8612814796873004 0.35095022356361616; 0.6884157053916526 0.481918985400
9987; -0.0435701516313082 0.8681552671858593; 0.7245715179296826 0.60524795
99696545]   
 [0.639416506953087 0.14149186863102034; 0.6045055288378626 0.5539149485415
189; -0.17604607024753233 0.8509308437858583; 1.0234926858206546 0.85811132
14950639]   
 [0.2045828895696477 -0.2525041978120318; 0.3189536789846515 0.498253464653
9205; -0.21859428329425323 0.9428260910353048; 1.3676916752649069 1.1293421
70650124]   
 [-0.5247863084275708 -0.8853558155271741; -0.22904535724602598 0.272081830
017429; -0.022586423462697486 1.2658051861544248; 1.6930321068954122 1.3473
66838750421]
 [-1.596241487638124 -1.7650129829900383; -1.0096364160682132 -0.0777405483
9070577; 0.6111815116390642 1.9710492033139315; 1.9013867799889466 1.410428
616184992]  
 [-3.088864730689933 -2.8884867482243854; -1.89939215748847 -0.377417741616
9476; 2.0255573851318265 3.3106758294172733; 1.8401389914264719 1.159055788
0944177]    
 [-4.741278898243673 -3.9335065059007217; -2.3837684701411317 -0.1628264011
8413036; 4.420117802509115 5.328489736198788; 1.271812824420278 0.396970463
19250806]   
 [-6.0300291310205925 -4.367029492312431; -1.74274192279828 1.1540427220849
248; 7.657029927596874 7.728063705591205; 0.00947891820170188 -0.9832882854
497728]     
 [-6.137700521949925 -3.2666820278834194; 1.3340496794841403 4.697017393957
745; 11.717909327425811 10.089825057925813; -2.4251766777177335 -3.31580501
33186405]   
 [-5.826048266206638 -2.666854071203864; 2.487015014508281 5.86393368970059
7; 12.557283203988641 10.417972995877923; -3.129665835379061 -3.94069849695
08274]
sol[3]
4×2 StaticArrays.SArray{Tuple{4,2},Float64,2,8} with indices SOneTo(4)×SOne
To(2):
  0.639417  0.141492
  0.604506  0.553915
 -0.176046  0.850931
  1.02349   0.858111

Conclusion

These are the basic controls in DifferentialEquations.jl. All equations are defined via a problem type, and the solve command is used with an algorithm choice (or the default) to get a solution. Every solution acts the same, like an array sol[i] with sol.t[i], and also like a continuous function sol(t) with a nice plot command plot(sol). The Common Solver Options can be used to control the solver for any equation type. Lastly, the types used in the numerical solving are determined by the input types, and this can be used to solve with arbitrary precision and add additional optimizations (this can be used to solve via GPUs for example!). While this was shown on ODEs, these techniques generalize to other types of equations as well.

Appendix

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using DiffEqTutorials
DiffEqTutorials.weave_file("introduction","01-ode_introduction.jmd")

Computer Information:

Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\accou\AppData\Local\atom\app-1.42.0\atom.exe"  -a
  JULIA_NUM_THREADS = 4

Package Information:

Status `~\.julia\dev\DiffEqTutorials\Project.toml`