sundials

Re-posted from: http://www.pkgupdate.com/?p=18

Introduction to Pkg Update

Welcome to Pkg Update. This new blog focuses on the Julia package ecosystem: how packages/organizations are changing, what new projects you should be aware of, and how Base issues are propagating through the package-sphere. I hope that this will be a nice summary of the vast changes that happen throughout the Julia ecosystem. I want this to be a resource which both package users and package developers could read to get a general understanding of the big movements in the community. The goal is to link together issues/PRs/discussions from the the vast number of repositories to give a coherent idea of what changes to anticipate and congratulate. As Base continues to slim (for those who don’t know, there’s a trend for Base to become slimmer for a leaner install, with more functionality provided by the Julia organizations such as JuliaMath), I hope this helps you stay up to date with where everything is going!

While I try my best to be active throughout the Julia-sphere, this is not something that can be done alone. Please feel free to contact me via email, Twitter, or Gitter to give a summary of what’s going on in your organization / area of expertise. Also, if anyone is willing to join the Pkg Update team, I would be happy to have you! I am looking for correspondents to represent different Julia organizations and to help summarize the changes.

That said, here’s a quick highlight of some recent updates both users and developers might need to know.

JuliaMath: Exciting Developments in LibM.jl

I have to give a shoutout to @musm for his incredible work in Libm.jl. This project started with a discussion on whether Julia should continue to use OpenLibm (Libm is the math library which supplies the approximations for common functions like “sin” and “exp”). From this discussion the ambitious idea of creating a pure-Julia Libm was put forward, and Libm.jl has continued with the momentum.

Sleef is a C library notable for it’s use of SIMD. Libm.jl is largely based on a port of that code to Julia, it also has some added safety and correctness for edge cases over the original code. While the ported code is based on the non-SIMD parts of the C source, it reliably takes advantage f SIMD hardware none-the-less — thanks to LLVM autovectorisation . The code is particularly vectorisation friendly, as it uses next-to-no branches or table lookups Most of the benchmarks for this implementation are within 2x of C programs calling Sleef. In practice, because of how lining and FFI works, this means that using the julia Sleef implementation is significantly faster than ccalling the C sleef library. Not to mention safer.

Other libraries (Musl and Cephes) now have partial ports to Julia as well. There is work being done to use the strategies of these libraries together to build an even better Libm which might be dubbed Amal.jl. This library, being a Julia implementation, will have much broader support for the type system, including native implementations for Float16 and quad-precision floats, leading to a much more robust mathematical foundation.

JuliaStats API Decision: Degrees of Freedom as “dof”

This issue and the associated PR went through rather quickly. It highlighted the fact that until now, degrees of freedom in the stats packages had a tendency to use “df”, which is a naming conflict with using “df” as an instance for a DataFrame (a very common usage). The decision was to change the naming for degrees of freedom to dof. Since this change is in StatsBase, you can expect that this convention will propagate throughout the statistics stack: “df” for dataframes and “dof” for degrees of freedom.

JuliaDiffEq: Upcoming Sundials API Update

Sundials.jl, the interface to the popular ODE suite, has recently undergone a large overhaul to help improve the wrapper and to stop memory leaks. This is being combined with major changes to the “simplified interface” and an upgrade in the Sundials version. Thus the upcoming tag will be a minor release with some breaking but exciting improvements, including the ability to use the Sundials ARKODE solvers.

JunoLab: The Plot Pane is Taking Off

With the Julia v0.5 release, Juno released a new version of their stack which included an updated plot pane. The developers were rallied and many plotting packages now support the plot pane. This includes support from Plots.jl, meaning that the plots will plot in the plot pane if the backend is compatible. While Plotly is not compatible with the plot pane, a recent PR has added plot pane compatibility to PlotlyJS. PyPlot and GR are compatible with the plot pane. Also, compatibility with Gadfly has been added. This means that, at least on the master version of the Juno stack and the plotting packages, the plot pane is all working! For regular users you may want to wait for this all to be released, but if you’re brave, it’s already ready for testing and I have been making good use of it. Note that in Plots.jl you can still use the `gui()` command to open up the non-plot pane window. This can be useful since some plotting backends (like PyPlot) are not interactive in the plot pane.

JuliaML Manifesto Released, and Major Progress Ensues

If you haven’t read Tom Brelof’s JuliaML manifesto, I would take a good look at it. He gives an introduction to JuliaML, the new machine learning based Julia organization, and shows how the Learn ecosystem is leading to an easy-to-use but highly extendable API. If you’re interested in machine learning, this is an organization to start following. Lots of progress is happening over here.

New Package Highlights

Lastly, I would like to finish off each post by highlighting some newer and promising projects. Since this is the first post and there’s lots to say, I picked a few.

ParallelDataTransfer.jl

(Disclaimer: this is one of my own libraries) ParallelDataTransfer.jl has hit the top of the 14-day change chart on Julia Pulse, rapidly picking up steam since it’s new release. The package is pretty simple: it offers a easy-to-use API for performing safe data transfers between processes. The README also shows how to use the library within type-stable functions, giving an easy way to construct highly parallel Julia algorithms.

Query.jl

Another package showing large movement on the Julia Pulse is Query.jl. This package allows you to query from Julia data sources, including any iterable, DataFrames, and Arrays. It even allows for queries through DataStreams sources, effectively connecting the queries to CSVs and SQLite databases. It’s not hard to see the promising future that has.

RNG.jl

A Google Summer of Code project which has caught my attention is RNG.jl. This package does exactly what you’d expect: it implements a bunch of random number generators. The benchmarks show that there are some promising RNGs which are both fast and suitable for parallel usage. I could see this library’s RNGs as becoming a staple for many Julia packages.

That’s the End of the First Post!

Let me know what you think of this new blog. I tried to cover as wide of a range as I could, but to get a better sense of the even larger community, I will need your help. Let me know if you’re willing to be an organization correspondent who can write one of many (multiple) dispatches, or if you don’t have the time, feel free to submit stories when possible.

Yesterday we had a look at Julia’s support for Calculus. The next logical step is to solve some differential equations. We’ll look at two packages today: Sundials and ODE.

The Sundials package is based on a library which implements a number of solvers for differential equations. First off you’ll need to install that library. In Ubuntu this is straightforward using the package manager. Alternatively you can download the source distribution.

sudo apt-get install libsundials-serial-dev

Next install the Julia package and load it.

julia> Pkg.add("Sundials")
julia> using Sundials

To demonstrate we’ll look at a standard “textbook” problem: a damped harmonic oscillator (mass on a spring with friction). This is a second order differential equation with general form

$\ddot{x} + a \dot{x} + b x = 0$

where $x$ is the displacement of the oscillator, while $a$ and $b$ characterise the damping coefficient and spring stiffness respectively. To solve this numerically we need to convert it into a system of first order equations:

$\begin{aligned} \dot{x} &= v \\ \dot{v} &= - a v - b x \end{aligned}$

We’ll write a function for those relationships and assign specific values to $a$ and $b$ .

julia> function oscillator(t, y, ydot)
           ydot[1] = y[2]
           ydot[2] = - 3 * y[1] - y[2] / 10
       end
oscillator (generic function with 2 methods)

Next the initial conditions and time steps for the solution.

julia> initial = [1.0, 0.0];                   # Initial conditions
julia> t = float([0:0.125:30]);                # Time steps

And finally use cvode() to integrate the system.

julia> xv = Sundials.cvode(oscillator, initial, t);
julia> xv[1:5,:]
5x2 Array{Float64,2}:
 1.0        0.0     
 0.97676   -0.369762
 0.908531  -0.717827
 0.799076  -1.02841 
 0.65381   -1.28741

The results for the first few time steps look reasonable: the displacement (left column) is decreasing and the velocity (right column) is becoming progressively more negative. To be sure that the solution has the correct form, have a look at the Gadfly plot below. The displacement (black) and velocity (blue) curves are 90° out of phase, as expected, and both gradually decay with time due to damping. Looks about right to me!

ODE

The ODE package provides a selection of solvers, all of which are implemented with a consistent interface (which differs a bit from Sundials).

julia> using ODE

Again we need to define a function to characterise our differential equations. The form of the function is a little different with the ODE package: rather than passing the derivative vector by reference, it’s simply returned as the result. I’ve consider the same problem as above, but to spice things up I added a sinusoidal driving force.

julia> function oscillator(t, y)
           [y[2]; - a * y[1] - y[2] / 10 + sin(t)]
       end
oscillator (generic function with 2 methods)

We’ll solve this with ode23(), which is a second order adaptive solver with third order error control. Because it’s adaptive we don’t need to explicitly specify all of the time steps, just the minimum and maximum.

julia> a = 1;                                  # Resonant
julia> T, xv = ode23(oscillator, initial, [0.; 40]);
julia> xv = hcat(xv...).';                     # Vector{Vector{Float}} -> Matrix{Float}

The results are plotted below. Driving the oscillator at the resonant frequency causes the amplitude of oscillation to grow with time as energy is transferred to the oscillating mass.

If we move the oscillator away from resonance the behavior becomes rather interesting.

julia> a = 3;                                  # Far from resonance

Now, because the oscillation and the driving force aren’t synchronised (and there’s a non-rational relationship between their frequencies) the displacement and velocity appear to change irregularly with time.

How about a double pendulum (a pendulum with a second pendulum suspended from its end)? This seemingly simple system exhibits a rich range of dynamics. It’s behaviour is sensitive to initial conditions, one of the characteristics of chaotic systems.

First we set up the first order equations of motion. The details of this system are explained in the video below.

julia> function pendulum(t, y)
           Y = [
           6 * (2 * y[3] - 3 * cos(y[1] - y[2]) * y[4]) / (16 - 9 * cos(y[1] - y[2])^2);
           6 * (8 * y[4] - 3 * cos(y[1] - y[2]) * y[3]) / (16 - 9 * cos(y[1] - y[2])^2)
           ]
           [
           Y[1];
           Y[2];
           - (Y[1] * Y[2] * sin(y[1] - y[2]) + 3 * sin(y[1])) / 2;
           - (sin(y[2]) - Y[1] * Y[2] * sin(y[1] - y[2])) / 2;
           ]
       end
pendulum (generic function with 1 method)

Define initial conditions and let it run…

julia> initial = [pi / 4, 0, 0, 0];            # Initial conditions -> deterministic behaviour
T, xv = ode23(pendulum, initial, [0.; 40]);

Below are two plots which show the results. The first is a time series showing the angular displacement of the first (black) and second (blue) mass. Next is a phase space plot which shows a different view of the same variables. It’s clear to see that there is a regular systematic relationship between them.

Next we’ll look at a different set of initial conditions. This time both masses are initially located above the primary vertex of the pendulum. This represents an initial configuration with much more potential energy.

julia> initial = [3/4 * pi, pi, 0, 0];         # Initial conditions -> chaotic behaviour

The same pair of plots now illustrate much more interesting behaviour. Note the larger range of angles, θ₂, achieved by the second bob. With these initial conditions the pendulum is sufficiently energetic for it to “flip over”. Look at the video below to get an idea of what this looks like with a real pendulum.

It’s been a while since I’ve played with any Physics problems. That was fun. The full code for today is available at github. Come back tomorrow when I’ll take a look at Optimisation in Julia.

The post #MonthOfJulia Day 21: Differential Equations appeared first on Exegetic Analytics.

juliabloggers.com

A Julia Language Blog Aggregator

Tag Archives: sundials

Pkg Update: Your Window to the Julia Package Ecosystem