Merge pull request #30 from ChrisRackauckas/diffeq

Switch to DifferentialEquations.jl

Author: dlfivefifty

.travis.yml

@@ -6,7 +6,6 @@ addons:
 os:
     - linux
 julia:
-    - 0.4
     - 0.5
     - nightly
 sudo: false
@@ -17,4 +16,3 @@ notifications:
 #    - julia -e 'Pkg.clone(pwd()); Pkg.build("RandomMatrices"); Pkg.test("RandomMatrices"; coverage=true)'
 after_success:
     - julia -e 'Pkg.add("Coverage"); cd(Pkg.dir("RandomMatrices")); using Coverage; Coveralls.submit(process_folder()); Codecov.submit(process_folder())'
-

README.md

@@ -3,7 +3,7 @@ RandomMatrices.jl
 
 Random matrix package for [Julia](http://julialang.org).
 
-[![RandomMatrices](http://pkg.julialang.org/badges/RandomMatrices_0.4.svg)](http://pkg.julialang.org/?pkg=RandomMatrices)
+[![RandomMatrices](http://pkg.julialang.org/badges/RandomMatrices_0.5.svg)](http://pkg.julialang.org/?pkg=RandomMatrices)
 [![Build Status](https://travis-ci.org/jiahao/RandomMatrices.jl.png?branch=master)](https://travis-ci.org/jiahao/RandomMatrices.jl)
 [![Coverage Status](https://coveralls.io/repos/jiahao/RandomMatrices.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/jiahao/RandomMatrices.jl?branch=master)
 [![codecov.io](https://codecov.io/github/jiahao/RandomMatrices.jl/coverage.svg?branch=master)](https://codecov.io/github/jiahao/RandomMatrices.jl?branch=master)

REQUIRE

@@ -1,6 +1,6 @@
-julia 0.4
+julia 0.5
 Combinatorics 0.2
 Distributions 0.8.10
-ODE 0.2.1
+OrdinaryDiffEq
 GSL 0.3.1
 Compat 0.9.2

demos/book/8/randomgrowth2.jl

@@ -1,6 +1,6 @@
 #randomgrowth2.jl
 using PyPlot
-using ODE
+using OrdinaryDiffEq
 
 function randomgrowth2()
     num_trials = 2000
@@ -32,16 +32,22 @@ function randomgrowth2()
         t0 = 4
         tn = -6
         dx = 0.005
-        deq = (t, y) -> [y[2]; t*y[1]+2*y[1]^3; y[4]; y[1]^2]
-        y0 = [airy(t0); airy(1,t0); 0; airy(t0)^2]      # boundary conditions
-        t, y = ode23(deq, y0, t0:-dx:tn)                # solve
-        F2 = Float64[exp(-y[i][3]) for i = 1:length(y)]                               # the distribution
-        f2 = gradient(F2, t)                            # the density
+        deq = function (t, y, dy)
+          dy[1] = y[2]
+          dy[2] = t*y[1]+2y[1]^3
+          dy[3] = y[4]
+          dy[4] = y[1]^2
+        end
+        y0 = big.([airy(t0); airy(1,t0); 0; airy(t0)^2]) # boundary conditions
+        prob = ODEProblem(deq,y0,(t0,tn))
+        sol = solve(prob,Vern8(), saveat=-dx, abstol=1e-12, reltol=1e-12) # solve
+        F2 = Float64[exp(-sol[i][3]) for i = 1:length(y)]  # the distribution
+        f2 = gradient(F2, t)                             # the density
 
         # add(p, Curve(t, f2, "color", "red", "linewidth", 3))
         # Winston.display(p)
         subplot(1,length(Ns),jj)
-        plot(t, f2, "r", linewidth = 3)
+        plot(sol.t, f2, "r", linewidth = 3)
         ylim(0, 0.6)
         println(mean(C))
     end

src/RandomMatrices.jl

@@ -3,7 +3,8 @@ module RandomMatrices
 using Combinatorics
 using Compat
 using GSL
-using ODE
+using OrdinaryDiffEq
+using DiffEqBase # This line is only needed on v0.5, and comes free from OrdinaryDiffEq on v0.6
 
 import Base: isinf, rand
 import Distributions: ContinuousUnivariateDistribution,

src/densities/TracyWidom.jl

@@ -27,7 +27,7 @@ immutable TracyWidom <: ContinuousUnivariateDistribution end
 Probability density function of the Tracy-Widom distribution
 
 Computes the Tracy-Widom distribution by directly solving the
-Painlevé II equation using the ode23 numerical integrator
+Painlevé II equation using the Vern8 numerical integrator
 
 # Arguments
 * `d::TracyWidom` or `Type{TracyWidom}`: an instance of `TracyWidom` or the type itself
@@ -38,39 +38,43 @@ function pdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0))
     t≤t0 && return 0.0
     t≥5  && return 0.0
 
-    ts, y = _solve_painleve_ii(t0, t)
+    sol = _solve_painleve_ii(t0, t)
 
-    ΔF2=exp(-y[end-1][3]) - exp(-y[end][3]) # the cumulative distribution
-    f2=ΔF2/(ts[end]-ts[end-1])              # the density at t
+    ΔF2=exp(-sol[end-1][3]) - exp(-sol[end][3]) # the cumulative distribution
+    f2=ΔF2/(sol.t[end]-sol.t[end-1])            # the density at t
 end
 pdf(d::Type{TracyWidom}, t::Real) = pdf(d(), t)
 
 """
 Cumulative density function of the Tracy-Widom distribution
 
 Computes the Tracy-Widom distribution by directly solving the
-Painlevé II equation using the ode23 numerical integrator
+Painlevé II equation using the Vern8 numerical integrator
 
 See `pdf(::TracyWidom)` for a description of the arguments.
 """
 function cdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0))
     t≤t0 && return 0.0
     t≥5  && return 1.0
-
-    ts, y = _solve_painleve_ii(t0, t)
-
-    F2=exp(-y[end][3])
+    sol = _solve_painleve_ii(t0, t)
+    F2=exp(-sol[end][3])
 end
 cdf(d::Type{TracyWidom}, t::Real) = cdf(d(), t)
 
 # An internal function which sets up the Painleve II differential equation and
-# runs it through the ode23 numerical integrator
+# runs it through the Vern8 numerical integrator
 function _solve_painleve_ii{S<:Real}(t0::S, t::S)
-    deq(t, y) = [y[2], t*y[1]+2y[1]^3, y[4], y[1]^2]
+    function deq(t, y, dy)
+        dy[1] = y[2]
+        dy[2] = t*y[1]+2y[1]^3
+        dy[3] = y[4]
+        dy[4] = y[1]^2
+    end
     a0 = airy(t0)
-    T = typeof(a0)
+    T = typeof(big(a0))
     y0=T[a0, airy(1, t0), 0, airy(t0)^2]    # Initial conditions
-    (ts, y)=ode23(deq, y0, [t0, t])         # Solve the ODE
+    prob = ODEProblem(deq,y0,(t0,t))
+    solve(prob, Vern8(), abstol=1e-12, reltol=1e-12)         # Solve the ODE
 end
 
 """

test/densities/TracyWidom.jl

@@ -7,7 +7,7 @@ using Base.Test
 @test cdf(TracyWidom, -10) == 0
 @test cdf(TracyWidom, 10) == 1
 
-if isdefined(:ODE) && isa(ODE, Module)
+if isdefined(:OrdinaryDiffEq) && isa(OrdinaryDiffEq, Module)
     t = rand()
     @test pdf(TracyWidom, t) > 0
     @test 0 < cdf(TracyWidom, t) < 1