The option fPIC and -shared are essential for Julia to be able to resolve the function in our library. Now we are almost ready to build our library. From the bash terminal we invoke:
Differential equations are ubiquitous throughout mathematics and the sciences. In fact, I myself have studied various forms of differential equations stemming from fields including biology, chemistry, economics, and climatology. What was interesting is that, although many different people are using differential equations for many different things, pretty much everyone wants the same thing: to quickly solve differential equations in their various forms, and make some pretty plots to describe what happened.
The goal of DifferentialEquations.jl is to do exactly that: to make it easy solve differential equations with the latest and greatest algorithms, and put out a pretty plot. The core idea behind DifferentialEquations.jl is that, while it is easy to describe a differential equation, they have such diverse behavior that experts have spent over a century compiling different ways to think about and handle differential equations. Most users will want to just brush past all of the talk about which algorithms simply ask: “I have this kind of differential equation. What does the solution look like?”
DifferentialEquations.jl’s User Interface
To answer that question, the user should just have to say what their problem is, tell the computer to solve it, and then tell the computer to plot it. In DifferentialEquations.jl, we use exactly those terms. Let’s look at an Ordinary Differential Equation (ODE): the linear ODE . It is described as the function
To use DifferentialEquations.jl, you first have to tell the computer what kind of problem you have, and what your data is for the problem. Recall the general ordinary differential equation is of the form
and initial condition , so in this case, we have an ODE with data and . DifferentialEquations.jl is designed as a software for a high-level language, Julia. There are many reasons for this choice, but the one large reason is its type system and multiple dispatch. For our example, we tell the machine what type of problem we have by building a DEProblem type. The code looks like this:
using DifferentialEquations
alpha = 0.5#Setting alpha to 1/2
f(y,t) = alpha*y
u0 = 1.5
prob = ODEProblem(f,u0)
where prob contains everything about our problem. You can then tell the computer to solve it and give you a plot by, well, solve and plot:
timespan = [0,1]# Solve from time = 0 to time = 1
sol = solve(prob,timespan)# Solves the ODE
plot(sol)# Plots the solution using Plots.jl
And that’s the key idea: the user should simply have to tell the program what the problem is, and the program should handle the details. That doesn’t mean that the user won’t have access to to all of the details. For example, we can control the solver and plotters in more detail, using something like
sol = solve(prob,alg=:RK4)# Unrolled and optimzed RK4
plot(sol,lw=3)# All of the Plots.jl attributes are available
However, in many cases a user may approach the problem for which they don’t necessarily know anything about the algorithms involved in approximating the problem, and so obfuscating the API with these names is simply confusing. One place where this occurs is solving stochastic differential equations (SDEs). These have been recently growing in popularity in many of the sciences (especially systems biology) due to their added realism and their necessity when modeling rare and random phenomena. In DifferentialEquations.jl, you can get started by simply knowing that an SDE problem is defined by the functions and in the form
g(u,t) = 0.3u
prob = SDEProblem(f,g,u0)
sol = solve(prob,timespan)
plot(sol)
If you wish to dig into the manual, you will see that the default solver that is used is a Rossler-SRI type of method and will (soon) have adaptivity which is complex enough to be a numerical analysis and scientific computing research project. And you can dig into the manual to find out how to switch to different solvers, but the key idea is that you don’t have to. Everything is simply about defining a problem, and asking for solutions and plots.
Julia was created to solve the many-language problem in scientific computing. Before people would have to write out the inner loops as C/Fortran, and bind it to a scripting language that was never designed with performance in mind. Julia has done extremely well as solving this problem via multiple-dispatch. Multiple dispatch is not just about ease of use, but it is also the core of what allows Julia to be fast . From a quote I am stealing from IRC: “Julia: come for the fast loops, stay for the type-system”.
In my view, the many-language problem always had an uglier cousin: the many-API problem. Every package has its own way of interacting with the user, and it becomes a burden to remember how all of them work. However, in Julia there seems to be a beautiful emergence of packages which solve the many-API problem via Julia’s multiple-dispatch and metaprogramming functionalities. Take for example Plots.jl. There are many different plotting packages in Julia. However, through Plots.jl, you can plot onto any “backend” (other plotting package) using just one API. You can mix and match plotting in PyPlot (matplotlib), GR, Plotly, and unicode. It’s all the same commands. Another example of this is JuMP. Its core idea is solver independence: you take your optimization problem, define the model in JuMP’s lingo, and then plug into many different solvers all by flipping a switch.
But why “Differential Equations”? Isn’t that broad?
Sure, there are packages for solving various types of differential equations, all specializing in one little part. But when I was beginning my PhD, quickly found that these packages were missing something. The different types of differential equations that we encounter are not different but many times embody the same problem: a PDE when discretized is a system of ODEs, the probability distribution of evolving SDEs is a PDE (a form of the Heat Equation), and all of the tools that they use to get good performance are the same. Indeed, many contemporary research questions can be boiled down to addressing the question: what happens if we change the type of differential equation? What happens if we add noise to our ODEs which describe population dispersal? What happens if we add to our model that RNA production is reading a delayed signal? Could we make this high-dimensional PDE computationally feasible via a Monte Carlo experiment combined with Feynman-Kac’s theorem?
Yet, our differential equations libraries are separate. Our PDEs are kept separate from our SDEs, while our delay equations hang out in their own world. Mixing and matching solvers requires learning complex APIs, usually large C/Fortran libraries with opaque function names. That is what DifferentialEquations.jl is looking to solve. I am building DifferentialEquations.jl as a hub for differential equations, the general sense of the term.
If you have defined an SDE problem, then via the Forward Kolmorogov equation there is a PDE associated to the SDE. In many cases like the Black-Scholes model, both the SDE and the PDE are canonical ways of looking at the same problem. The solver should translate between them, and the solver should handle both types of differential equations. With one API and the functionality for these contained within the same package, no longer are they separate entities to handle computationally.
Where are we currently?
DifferentialEquations.jl is still very young. Indeed, the project only started a few months ago, and during that time period I was taking 6 courses. However, the package already has a strong core, including
You may have been thinking, “but I am a numerical analyst. How could this program help me?”. DifferentialEquations.jl has a bunch of functionalities for quickly designing and testing algorithms. All of the DEProblems allow for one to give them the analytical solution, and the solvers will then automatically calculate the errors. Thus by using some simple macros, one can define new algorithms in just a few lines of code, test the convergence, benchmark times, and have the algorithm available as an `alg` option in no time (note: all of the ODE solvers were written in one morning!). Thus it is easy to define the loop, and the rest of the functionality will come by default. It’s both a great way to test algorithms, and share algorithms. Contributing will both help you and DifferentialEquations.jl!.
Where are we going?
I have big plans for DifferentialEquations.jl. For example:
I will be rolling out an efficiency testing suite so that one could just specify the algorithms you’d like to solve a problem, and along what convergence axis (i.e. choose a few s, or along changing tolerances), and it will output comparisons of the computational efficiencies and make some plots. It will be similar in functionality to the ConvergenceSimulation suite.
Finite difference methods for Heat and Poisson equation. These are long overdue for the research I do.
Changing the tableaus in ODEs and SDEs to StaticArrays so they are stack allocated. This has already been tested and works well on v0.5.
Higher-order methods for parabolic SPDEs (a research project with promising results!).
Blazing fast adaptivity for SDEs. (Once the paper I have submitted for it is published, it will be available. It’s already implemented!)
High-stability high order methods for SDEs (another research project).
Parallel methods. I have already implemented parallel (Xeon Phi) solvers and described them in previous blog posts. They simply need to be integrated into DifferentialEquations.jl. I would like to have native GPU solvers as well.
Delay and algebraic differential equations.
Wrapping more popular solvers. I’d like to add Sundials, LSODE, and PetsC to the list.
A web interface via Escher.jl to define DEProblems and get the solution plots. I am looking to have this hosted as an XSEDE Gateway.
If you’d like to influence where this project is going, please file an issue on the Github repository. I am always open for suggestions.
I hope this gives you a good idea on what my plans are for DifferentialEquations.jl. Check out the documentation and give the package a whirl!
So far in this blog, all my plots have either used PyPlot, or the experimental GLVisualize. BUT, NO LONGER!
While writing up my research notes, I’ve been trying to improve my PyPlot graphs. I was hoping for them to look nicer; Julia’s PyPlot was saying no. We had a falling out over background colors and resolution. So, I have decided to move on to a plotting partner more willing to cope with my demands.
When I first started learning Julia, Gadfly seemed too unfamiliar, and control too tricky. I didn’t give it much further thought. After some digging, Gadfly has changed my mind, and I hope today I can change yours too. Right now, like Julia, Gadfly doesn’t have the maturity or documentation of PyPlot, but once it does, you better watch out.
Gadfly preferentially produces plots in SVG (scalable vector graphics) form, instead of the rasterized form of PNG or JPG favored by PyPlot. As displayed in the figure left, rasterized images are composed of colored squares, so don’t look the best when highly magnifified. SVG, which is based on XML, stores the position of vector elements and scales to high resolution without artifacts. Tools such as Inkscape also work in SVG, and can edit the individual pieces of a graph since they haven’t been rasterized together.
Let’s plot some graphs!
By Yug, modifications by 3247 (Unknown) [CC BY-SA 2.5 (http://creativecommons.org/licenses/by-sa/2.5)], via Wikimedia Commons
</div>
# Time to load our Packages#Pkg.update()#Pkg.add("Gadfly")usingGadfly#Pkg.add("Compose")usingCompose
# Some test datax=collect(1:.4:10);y1=sin(x);y2=cos(x);
# A simple plotplot(x=x,y=y1)
Unlike PyPlot, we can give plot functions as well as points. After receiving a function, a max value, and a min value, Gadfly figures everything else out for itself.
For the function, we can pass in either an inbuilt function, a user defined function, or an anonymous function.
The brackets [ ] allow us to group multiple functions together in the plot.
But that’s just one set of points…
How do we plot more than one set of data? That’s where layers come in. Instead of using plot, we use layer, and set it as a variable. We can then plot those layers.
Now, what if we want to plot something other than lines?
Gadfly allows a wide variety of other ways to plot our data. We control this with the Geom (Geometry) object.
Disclaimer: I’ve learned how to do the next two things by reading the code, not the documentation.
It seems like they are newer and less developed features than everything else I’m discussing here. The syntax seems less polished than in other areas, so I believe it’s still under construction.
# or Compose.mm for smaller sizes# These ratios and numbers changed around how ever you likedash=6*Compose.cmdot=1.5*Compose.cmgap=1*Compose.cml1=layer(x=x,y=zeros(x),Geom.line,Theme(line_style=[dot]))l2=layer(x=x,y=ones(x),Geom.line,Theme(line_style=[dash]))l3=layer(x=x,y=2*ones(x),Geom.line,Theme(line_style=[gap]))l4=layer(x=x,y=3*ones(x),Geom.line,Theme(line_style=[dot,dash,dot]))l5=layer(x=x,y=4*ones(x),Geom.line,Theme(line_style=[gap,dash]))plot(l1,l2,l3,l4,l5,Coord.Cartesian(ymin=-.5,ymax=4.5),Theme(grid_line_width=0pt))
Guide: Labeling Axes
Where Geom alters how we plot, Guide alters the labeling. Guide ties in with Compose.jl through the Guide.annotate command, but that will take a tutorial in itself.
Here’s something we can do with a combination of Guide and Scale. Using Guide, we can set where we specifically want our xticks to be, namely at multiples of $\pi$. But then, the labeling would write out some irrational number, making the plot look horrible. So, we create a function that takes in a number and outputs a string label for Scale to use.
function labelname(x::Real)n=round(Int,x/π)#nearest integer*piifn==0return"0"elseifn==1return"π"endreturn("$n π")endxmarks=[0,π,2π,3π]ymarks=[-1,0,1]plot(x=x,y=y2,Guide.yticks(ticks=ymarks),Guide.xticks(ticks=xmarks),Scale.x_continuous(labels=labelname))
Some other cool things we can do with Scale:
* Automatically transform the axis according to log, log10,log2, asinh, or sqrt.
* Write numbers in :plain, :scientific, :engineering, or :auto
I mostly chose Gadfly because of the control I could have with the Theme command. http://dcjones.github.io/Gadfly.jl/themes.html has a much more exhaustive list than what I will be demonstrating with here.
# Solarized Colors that I like working with# http://ethanschoonover.com/solarizedusingColorsbase03=parse(Colorant,"#002b36");base02=parse(Colorant,"#073642");base01=parse(Colorant,"#586e75");base00=parse(Colorant,"#657b83");base0=parse(Colorant,"#839496");base1=parse(Colorant,"#839496");base2=parse(Colorant,"#eee8d5");base3=parse(Colorant,"#fdf6e3");yellow=parse(Colorant,"#b58900");orange=parse(Colorant,"#cb4b16");red=parse(Colorant,"#dc322f");magenta=parse(Colorant,"#d33682");violet=parse(Colorant,"#6c71c4");blue=parse(Colorant,"#268bd2");cyan=parse(Colorant,"#3aa198");green=parse(Colorant,"#859900");
plot(x=x,y=y1,Theme(highlight_width=0pt,# lose the white ring around the pointsdefault_point_size=3pt,# size of the dotsdefault_color=magenta,# color of the dotsbackground_color=base03,# ... background color ...grid_color=base2,# the linesminor_label_color=base2,# numbers on axes colorkey_label_color=base2,# color key numbers colorkey_title_color=cyan,# label of color key colormajor_label_color=cyan,# axes label and title colormajor_label_font_size=20pt,# font size for axes label and titleminor_label_font_size=15pt,#font size for numbers on axespanel_opacity=1#setting background to opaque);)
Coord: Setting the boundries
The documentation states that you can change the range of the axes using Scale, but I’ve found it better to use this format to set my min and max values.
So far, we have covered seperately Guide, Themes, and partially Coord and Scale. Individually, each aspect doesn’t add too much clunkiness to the code. However, if we start to add everything together, then the plot function would look quite nasty.
Luckily, just like layers for data points, we can put our modifiers into variables. Then we can simply call the variables in our plot function.
This also helps for when we want to use one theme for a series of graphs. We can define the theme variable up top, and then include it in every graph there after, never having to declare it again. This helps me to keep my graphs uniform in style.
function labelname(x::Real)n=round(Int,x/π)ifn==0return"0"elseifn==1return"π"elsereturn("$n π")endreturn("$n π")endxticks=[0,π,2π,3π]yticks=[-1,0,1]data=layer(x=x,y=y1,ymin=y1-.1,ymax=y1+.1,Geom.line,Geom.ribbon)f=layer(cos,0,3π)yt=Guide.yticks(ticks=yticks)xt=Guide.xticks(ticks=xticks)xs=Scale.x_continuous(labels=labelname)t=Theme(highlight_width=0pt,default_point_size=3pt,default_color=blue,background_color=base03,grid_color=base2,minor_label_color=base2,key_label_color=base2,key_title_color=cyan,major_label_color=cyan,major_label_font_size=20pt,minor_label_font_size=15pt,panel_opacity=1)xl=Guide.xlabel("time")yl=Guide.ylabel("y axis")GT=Guide.title("Important title")plot(data,f,yt,xt,xs,t,xl,yl,GT)
Although we still have to call quite a few variables, this is a much simpler way of doing things.
Saving the Figure
Julia naturally saves to SVG (or SVGJS).
We have to specify the x and y dimensions of the plot, but since these images rescale so easily, we can choose some reasonable numbers.
