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.9730384275622996
 2.664456142481451

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.9730384275622996
 2.664456142481451

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.66648203038311
 3.664456142481451

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.9730384275622996
 2.664456142481451

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.554261048055312

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.102962785129292
 1.2165269512238264
 1.341783821227542
 1.4799379510586077
 1.632316207054161
 1.8003833264983584
 1.9857565541588758
 2.1902158127997695
 2.415725742084496
 2.664456142481451

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.9857565541588758
 2.415725742084496

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.9730384275622996
 2.664456142481451

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.664456142481451

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.664456142481451

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.0434491434754793
 2.641825616057761
 2.6644526430553817

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.090571742234628
 2.486102171447025
 2.6644562434913377

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: 10-element Array{Float64,1}:
 0.0
 0.05065844961708149
 0.12892424499087696
 0.2202665346951705
 0.33536905471453515
 0.45849793339130746
 0.6034831469120607
 0.7571559089306192
 0.9286657004231316
 1.0
u: 10-element Array{Array{Float64,2},1}:
 [0.38463437652876764 0.3365394358001754; 0.95367647369539 0.90977595277039
12; 0.4933912443502877 0.2725101654104203; 0.9517662367605697 0.11604197538
996264]
 [0.13320745263393463 0.3191628244842401; 0.8549693932660031 0.912536712052
0655; 0.4928856399687078 0.21271418893558364; 1.1347220680936703 0.15174973
68095521]
 [-0.3655044182392014 0.274059779139162; 0.5837099646311803 0.8778237974396
547; 0.6240289604081771 0.1330054287582697; 1.3585881373075774 0.1957014562
3446372]
 [-1.0900211249354295 0.1982869571314106; 0.14963513954964586 0.78869942082
94823; 1.0176918742754268 0.06556094183410746; 1.5041941575478772 0.2276541
6557284213]
 [-2.141907187516936 0.07941000409647822; -0.41785950728539584 0.6279099916
099771; 1.9313235255215078 0.027942899162192864; 1.4615974058994936 0.23536
809043318377]
 [-3.2768398098633433 -0.05594528651475372; -0.7812392978462199 0.436964469
2408842; 3.42979980463502 0.0496819173052162; 1.0744392775521392 0.20166553
467571022]
 [-4.291879866096152 -0.18998313793950897; -0.4377885447102866 0.2498866714
4071008; 5.743546563885415 0.14645875431949915; 0.09484743526138606 0.10804
830794775776]
 [-4.447324728841423 -0.2524569296692268; 1.3817233458419276 0.166354759676
60782; 8.388009798344632 0.2934187694456827; -1.5877443480758702 -0.0455685
3527350534]
 [-2.6856298057401657 -0.1624660959127634; 5.6709708257812865 0.27178161889
36869; 10.494669977060056 0.42139185264180334; -4.121425126129509 -0.258028
71833301305]
 [-1.1477516885708174 -0.062387615575211317; 8.17413480570777 0.38035829861
31633; 10.718557987022905 0.43274694875320296; -5.290842001934896 -0.349054
2998821476]

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.365504  0.27406
  0.58371   0.877824
  0.624029  0.133005
  1.35859   0.195701

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.384634  0.336539
 0.953676  0.909776
 0.493391  0.27251
 0.951766  0.116042

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: 5-element Array{Float64,1}:
 0.0
 0.19070503676179112
 0.4856523153591445
 0.8019021954559332
 1.0
u: 5-element Array{Array{BigFloat,2},1}:
 [0.3846343765287676408348715995089150965213775634765625 0.3365394358001754
238074454406159929931163787841796875; 0.95367647369538999235771825624397024
51229095458984375 0.90977595277039124965767769026570022106170654296875; 0.4
933912443502876943313140145619399845600128173828125 0.272510165410420279741
2881409400142729282379150390625; 0.9517662367605697060213287841179408133029
937744140625 0.11604197538996263716626344830729067325592041015625]
 [-0.8413141448946438464854200006224046557271134002128273977497375649710760
146303335 0.225115946417345641108203738065506990283035143541402544071612504
1178357579627907; 0.2983934653892308448851591887145020577821645180731939732
798287766022537113105765 0.822338578014582963730232151119317339003077097205
0389887974037912623200875271187; 0.8595861341673985635459003973116980641004
58442998669611869656295993803680928434 0.0839947462932442013741223373114517
9676304048175674375038692263909912292159854536; 1.4725869303236413502141549
83113092066857603617934887626247419036627068380840933 0.2197241087609361405
6997980818413493559961866839068810782958803548548682928695]
 [-3.5056462528545846058441423923326503886772387531077154966736064351766321
54944261 -0.084358133144244691300541350448192767513651585460994401259968861
77943322850398993; -0.79493765518031331455452601041382473598230260810489201
43579017783465281654159142 0.3968709158942837879409717701557494882133281293
130540231546992471877343588187838; 3.82560191316957853794568064799022432838
6747805315186610957782931562880447299358 0.06261201320113883738257842626144
93442695586472530594341404198240080097622066548; 0.935731961602460588705342
682373790191440147493320915573316023570417936846163878 0.188434429324558619
8122453414351697785018639346573316828021023041007501927339766]
 [-4.2197452336612217997217348665326508876419431016624007174140771887483192
76141243 -0.247654305990673897792364034877336387644193481432213534570862342
7783746689644428; 2.2598413470265915799026316286121460034343355660481614322
94891839011465928156635 0.1722990551198694989166638204189549965735630226076
652639607644534809934568272778; 9.08142317380597448358782610370147777225939
5397302976982078091594101492534768479 0.33514699620535340564036404277211145
37551218395670101375701854453584367074603045; -2.19273072389578513496240355
9581985020992175013626845744558150005543535712908488 -0.0982525156740925663
2586861519445680147932379077905944940285751065325873136563135]
 [-1.1477353214836802025571388330521720415737497614680557697242172358296974
90012226 -0.062386622795402500138936845062374804749730881829237525837944184
91099708538729871; 8.174174560295296718552076839071227804555870265831988390
156469810640243062081908 0.380359883974219886908715403677224665172803864370
1478681610840763735504735805074; 10.718572377696180607107634922956243258440
67453465441143710735830361224161528317 0.4327470606444676352017784320570234
875995526345834662868897024262073516028952793; -5.2908615329812109654432443
15945928586878501042216582713558989971465271793996808 -0.349055606103435442
0692200940818770604278331987750973184454922050420485595339426]
sol[1,3]
-3.505646252854584605844142392332650388677238753107715496673606435176632154
944261

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: 5-element Array{BigFloat,1}:
 0.0
 0.190705036761791123761979678319011593271217634692862377355564380641419549
1509929
 0.485652315359144500997473437435559286036998845479595580077598767278892407
2433831
 0.801902195455933254952655589517056579554309960231940038944772964506178528
1607648
 1.0
u: 5-element Array{Array{BigFloat,2},1}:
 [0.3846343765287676408348715995089150965213775634765625 0.3365394358001754
238074454406159929931163787841796875; 0.95367647369538999235771825624397024
51229095458984375 0.90977595277039124965767769026570022106170654296875; 0.4
933912443502876943313140145619399845600128173828125 0.272510165410420279741
2881409400142729282379150390625; 0.9517662367605697060213287841179408133029
937744140625 0.11604197538996263716626344830729067325592041015625]
 [-0.8413141448946438394289549426204327688468748210693597601341946160991893
933531594 0.225115946417345641859504048703528530244615484110736890447943115
626433873575101; 0.29839346538923084913530310648123256101651982829022392611
06747579301899978461498 0.8223385780145829646483074491670013533820940420891
392304011294308131792464652417; 0.85958613416739855938487669511421993730682
39019555091556737106262774668498676812 0.0839947462932442019596257597074876
6518047977240339776820785445155661273048297443; 1.4725869303236413490671482
95902459477751610238940088828287274202947321321058838 0.2197241087609361403
050139723243169872173142082768603424612937222459640527599299]
 [-3.5056462528545846711983305447101614720209067279619273387274444163129089
41147019 -0.084358133144244699497699992325594017329834062713710971634455530
91147812526017127; -0.79493765518031331405542163863147201872949937818027445
25022819391417938072679646 0.3968709158942837763939825317285620355056716111
658046192786168455326975597132107; 3.82560191316957865739242458654913376631
9094690522621015675084222260774073991841 0.06261201320113884158178689805210
702590450643350866799729579275750412871134401055; 0.93573196160246054494621
07364690305333527266000139681788251984899089718343876132 0.1884344293245586
156199180114035446946984710096114292439494295500966943147163976]
 [-4.2197452336612215635002086152297664800693348035626437780046957335605586
88091396 -0.247654305990673889259400539606246997771581500316147824183661337
4112775654100622; 2.2598413470265923331743041050478058092816374099716175189
05572391919167617388353 0.1722990551198695094294392362009935940431221093169
58486516601559386742874460678; 9.081423173805974998008677173025578602327863
223280758428480859043913400069206162 0.335146996205353436897559161327493923
8638567986310656890500111301895101679672078; -2.192730723895785626528117078
201929006853113297892289185702214668364928757972765 -0.09825251567409260861
559154967720811337157239243997321510423081021096705758806748]
 [-1.1477353214836802025571001641656068352507943452523479639133071423878665
39426396 -0.062386622795402500138935673281468942269581979795828624120814136
36505305755715494; 8.174174560295296718552037246925171090802501913212094753
083429952790930156581072 0.380359883974219886908710433749077779164302977325
3702111165378748036978107212245; 10.718572377696180607107516939616896463071
02431425663353757038486620605184074456 0.4327470606444676352017703616803117
227836778433138859357524055119484757259852219; -5.2908615329812109654432088
09709306876235067140433262050038794191804493100698004 -0.349055606103435442
0692173301891463825658509462390173904559580222993975422700265]

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: 10-element Array{Float64,1}:
 0.0
 0.048006104540199344
 0.12508535396737538
 0.2128933453419654
 0.32374584194942524
 0.4446150526079904
 0.5833187308680358
 0.7299931730760287
 0.8834638751632482
 1.0
u: 10-element Array{StaticArrays.SArray{Tuple{4,2},Float64,2,8},1}:
 [0.6470268048115251 0.42689239938470225; 0.4468712596313109 0.930792622748
21; 0.8995497299530391 0.4955864508508838; 0.19897099052479028 0.7955219510
043468]
 [0.602392813583938 0.2316510143444062; 0.6326231722354123 0.87331931673791
55; 0.7936623605638815 0.4738199645556346; 0.4236254028720524 0.96137318774
05786]
 [0.41258139268512584 -0.1804019440562041; 0.769282141598299 0.668409010558
5981; 0.6794213884975377 0.5458882638846405; 0.7637205865686998 1.180165833
722606]
 [0.02194072683165854 -0.7782737762342137; 0.7054948276902163 0.31552386775
256436; 0.6800307037504818 0.8220189345947837; 1.0966484744079248 1.3372547
190895772]
 [-0.7101938625822555 -1.6694784653940866; 0.3716562266183081 -0.1881217526
428317; 0.9628967051367382 1.5119594167855097; 1.3868615058813365 1.3554355
28105561]
 [-1.7332757050882404 -2.6890904446257506; -0.13852866977255995 -0.59266048
01359094; 1.721984075196057 2.718169853049194; 1.4690378922944827 1.0921911
459488223]
 [-3.0034987304271112 -3.666588030802566; -0.5701335090616335 -0.5102704207
897966; 3.2257135247869404 4.604820390072593; 1.1673356535511101 0.36304997
29154003]
 [-4.10323291282484 -4.066217941561766; -0.34476360963855546 0.661994078047
2835; 5.448360327248404 6.881372397693559; 0.2934308635954497 -0.9404065360
691877]
 [-4.437996394845046 -3.226698038570903; 1.2302814174986403 3.4921094452365
8; 8.06195674098168 8.917139630712603; -1.2783527190357418 -2.8386164247715
86]
 [-3.726546412830267 -1.367370783576439; 3.63590816493634 6.832838070121893
; 9.765012397461916 9.606747074260191; -2.8823394142024394 -4.5229393664650
89]
sol[3]
4×2 StaticArrays.SArray{Tuple{4,2},Float64,2,8} with indices SOneTo(4)×SOne
To(2):
 0.412581  -0.180402
 0.769282   0.668409
 0.679421   0.545888
 0.763721   1.18017

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 SciMLTutorials.jl repository, found at: https://github.com/SciML/SciMLTutorials.jl. For more information on doing scientific machine learning (SciML) with open source software, check out https://sciml.ai/.

To locally run this tutorial, do the following commands:

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

Computer Information:

Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_LOAD_PATH = /builds/JuliaGPU/DiffEqTutorials.jl:
  JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqTutorials.jl/.julia
  JULIA_CUDA_MEMORY_LIMIT = 2147483648
  JULIA_NUM_THREADS = 8

Package Information:

Status `/builds/JuliaGPU/DiffEqTutorials.jl/tutorials/introduction/Project.toml`
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.15.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 5.4.0
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.5.8
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.12.4
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 4.2.5
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra