From c2de0291a864799273ff6f59864812e613cd8f13 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 16:47:13 +0200 Subject: [PATCH 01/14] Component-wise RelTol and AbsTol for Explicit Solvers similar to Solving Ordinary Differential Equations I by Hairer and Wanner. Component-wise AbsTol for Implicit Solvers similar to MATLAB ode23s. --- .ipynb_checkpoints/Untitled-checkpoint.ipynb | 487 +++++++++++++++ .../Untitled-checkpoint.ipynb | 582 ++++++++++++++++++ src/ODE.jl | 32 +- src/Untitled.ipynb | 582 ++++++++++++++++++ src/runge_kutta.jl | 9 +- test/runtests.jl | 2 +- 6 files changed, 1681 insertions(+), 13 deletions(-) create mode 100644 .ipynb_checkpoints/Untitled-checkpoint.ipynb create mode 100644 src/.ipynb_checkpoints/Untitled-checkpoint.ipynb create mode 100644 src/Untitled.ipynb diff --git a/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/.ipynb_checkpoints/Untitled-checkpoint.ipynb new file mode 100644 index 000000000..97d2e1d78 --- /dev/null +++ b/.ipynb_checkpoints/Untitled-checkpoint.ipynb @@ -0,0 +1,487 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "LoadError", + "evalue": "LoadError: could not open file /mnt/526E44596E4437CD/GitRepo/ODE.jl/runge_kutta.jl\nwhile loading In[1], in expression starting on line 149", + "output_type": "error", + "traceback": [ + "LoadError: could not open file /mnt/526E44596E4437CD/GitRepo/ODE.jl/runge_kutta.jl\nwhile loading In[1], in expression starting on line 149", + "", + " in include at ./boot.jl:261", + " in include_from_node1 at ./loading.jl:304" + ] + } + ], + "source": [ + "#isdefined(Base, :__precompile__) && __precompile__()\n", + "# Ordinary Differential Equation Solvers\n", + "\n", + "#module ODE\n", + "\n", + "using Polynomials\n", + "using Compat\n", + "\n", + "## minimal function export list\n", + "# adaptive non-stiff:\n", + "export ode23, ode45, ode78\n", + "# non-adaptive non-stiff:\n", + "export ode4, ode4ms\n", + "# adaptive stiff:\n", + "export ode23s\n", + "# non-adaptive stiff:\n", + "export ode4s\n", + "\n", + "## complete function export list: see runtests.jl\n", + "\n", + "###############################################################################\n", + "## Coefficient Tableaus\n", + "###############################################################################\n", + "\n", + "# Butcher Tableaus, or more generally coefficient tables\n", + "# see Hairer & Wanner 1992, p. 134, 166\n", + "\n", + "abstract Tableau{Name, S, T<:Real}\n", + "# Name is the name of the tableau/method (a symbol)\n", + "# S is the number of stages (an int)\n", + "# T is the type of the coefficients\n", + "#\n", + "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", + "#\n", + "# For all types of tableaus it assumes fields:\n", + "# order::(Int...) # order of the method(s)\n", + "#\n", + "# For Runge-Kutta methods it assumes fields:\n", + "# a::Matrix{T} # SxS matrix\n", + "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", + "# c::Vector{T} # S\n", + "#\n", + "# For a tableau:\n", + "# c1 | a_11 .... a_1s\n", + "# . | a_21 . .\n", + "# . | a_31 . .\n", + "# . | .... . .\n", + "# c_s | a_s1 ....... a_ss\n", + "# -----+--------------------\n", + "# | b_1 ... b_s this is the one used for stepping\n", + "# | b'_1 ... b'_s this is the one used for error-checking\n", + "\n", + "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", + "order(b::Tableau) = b.order\n", + "# Subtypes need to define a convert method to convert to a different\n", + "# eltype with signature:\n", + "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", + "\n", + "###############################################################################\n", + "## HELPER FUNCTIONS\n", + "###############################################################################\n", + "\n", + "# estimator for initial step based on book\n", + "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", + "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", + " # Returns first step, direction of integration and F evaluated at t0\n", + " tdir = sign(tend-t0)\n", + " tdir==0 && error(\"Zero time span\")\n", + " tau = max(reltol*norm(x0, Inf), abstol)\n", + " d0 = norm(x0, Inf)/tau\n", + " f0 = F(t0, x0)\n", + " d1 = norm(f0, Inf)/tau\n", + " if d0 < 1e-5 || d1 < 1e-5\n", + " h0 = 1e-6\n", + " else\n", + " h0 = 0.01*(d0/d1)\n", + " end\n", + " # perform Euler step\n", + " x1 = x0 + tdir*h0*f0\n", + " f1 = F(t0 + tdir*h0, x1)\n", + " # estimate second derivative\n", + " d2 = norm(f1 - f0, Inf)/(tau*h0)\n", + " if max(d1, d2) <= 1e-15\n", + " h1 = max(1e-6, 1e-3*h0)\n", + " else\n", + " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", + " h1 = 10.^pow\n", + " end\n", + " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", + "end\n", + "\n", + "# isoutofdomain takes the state and returns true if state is outside\n", + "# of the allowed domain. Used in adaptive step-control.\n", + "isoutofdomain(x) = isnan(x)\n", + "\n", + "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", + " # There are a few types involved in a call to a ODE solver which\n", + " # somehow need to be consistent:\n", + " #\n", + " # Et = eltype(tspan)\n", + " # Ey = eltype(y0)\n", + " # Ef = eltype(Tf)\n", + " #\n", + " # There are also the types of the containers, but they are not\n", + " # needed as `similar` is used to make containers.\n", + " # Tt = typeof(tspan)\n", + " # Ty = typeof(y0) # note, this can be a scalar\n", + " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", + " #\n", + " # Returns\n", + " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", + " # the moment a AbstractFloat\n", + " # - Eyf: suitable eltype of y and f(t,y)\n", + " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", + " # - Ty: container type of y0\n", + " # - btab: tableau with entries converted to Et\n", + "\n", + " # Needed interface:\n", + " # On components: /, -\n", + " # On container: eltype, promote_type\n", + " # On time container: eltype\n", + "\n", + " Ty = typeof(y0)\n", + " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", + "\n", + " Et = eltype(tspan)\n", + " @assert Et<:Real\n", + " if !(Et<:AbstractFloat)\n", + " Et = promote_type(Et, Float64)\n", + " end\n", + "\n", + " # if all are Floats, make them the same\n", + " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", + " Et = promote_type(Et, Eyf)\n", + " Eyf = Et\n", + " end\n", + "\n", + " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", + " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", + "\n", + " btab_ = convert(Et, btab)\n", + " return Et, Eyf, Ty, btab_\n", + "end\n", + "\n", + "###############################################################################\n", + "## NON-STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "include(\"runge_kutta.jl\")\n", + "\n", + "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", + "# with Adams-Bashforth-Moulton coefficients\n", + "function ode_ms(F, x0, tspan, order::Integer)\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " if 1 <= order <= 4\n", + " b = ms_coefficients4\n", + " else\n", + " b = zeros(order, order)\n", + " b[1:4, 1:4] = ms_coefficients4\n", + " for s = 5:order\n", + " for j = 0:(s - 1)\n", + " # Assign in correct order for multiplication below\n", + " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", + " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", + " b[s, s - j] = ((-1)^j / factorial(j)\n", + " / factorial(s - 1 - j) * polyval(p_int, 1))\n", + " end\n", + " end\n", + " end\n", + "\n", + " # TODO: use a better data structure here (should be an order-element circ buffer)\n", + " xdot = similar(x)\n", + " for i = 1:length(tspan)-1\n", + " # Need to run the first several steps at reduced order\n", + " steporder = min(i, order)\n", + " xdot[i] = F(tspan[i], x[i])\n", + "\n", + " x[i+1] = x[i]\n", + " for j = 1:steporder\n", + " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", + " end\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "# Use order 4 by default\n", + "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", + "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", + "\n", + "###############################################################################\n", + "## STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "# Crude forward finite differences estimator of Jacobian as fallback\n", + "\n", + "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", + "function fdjacobian(F, x::Number, t)\n", + " ftx = F(t, x)\n", + "\n", + " # The 100 below is heuristic\n", + " dx = (x .+ (x==0))./100\n", + " dFdx = (F(t,x+dx)-ftx)./dx\n", + "\n", + " return dFdx\n", + "end\n", + "\n", + "function fdjacobian(F, x, t)\n", + " ftx = F(t, x)\n", + " lx = max(length(x),1)\n", + " dFdx = zeros(eltype(x), lx, lx)\n", + " for j = 1:lx\n", + " # The 100 below is heuristic\n", + " dx = zeros(eltype(x), lx)\n", + " dx[j] = (x[j] .+ (x[j]==0))./100\n", + " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", + " end\n", + " return dFdx\n", + "end\n", + "\n", + "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", + "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", + "#\n", + "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", + "#\n", + "# supports keywords: points = :all | :specified (using dense output)\n", + "# jacobian = G(t,y)::Function | nothing (FD)\n", + "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", + " jacobian=nothing,\n", + " points=:all,\n", + " norm=Base.norm,\n", + " minstep=abs(tspan[end] - tspan[1])/1e18,\n", + " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", + " initstep=0.)\n", + "\n", + "\n", + " # select method for computing the Jacobian\n", + " if typeof(jacobian) == Function\n", + " jac = jacobian\n", + " else\n", + " # fallback finite-difference\n", + " jac = (t, y)->fdjacobian(F, y, t)\n", + " end\n", + "\n", + " # constants\n", + " const d = 1/(2 + sqrt(2))\n", + " const e32 = 6 + sqrt(2)\n", + "\n", + "\n", + " # initialization\n", + " t = tspan[1]\n", + "\n", + " tfinal = tspan[end]\n", + "\n", + " h = initstep\n", + " if h == 0.\n", + " # initial guess at a step size\n", + " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", + " else\n", + " tdir = sign(tfinal - t)\n", + " F0 = F(t,y0)\n", + " end\n", + " h = tdir * min(abs(h), maxstep)\n", + "\n", + " y = y0\n", + " tout = Array(typeof(t), 1)\n", + " tout[1] = t # first output time\n", + " yout = Array(typeof(y0), 1)\n", + " yout[1] = deepcopy(y) # first output solution\n", + "\n", + "\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + "\n", + " while abs(t - tfinal) > 0 && minstep < abs(h)\n", + " if abs(t-tfinal) < abs(h)\n", + " h = tfinal - t\n", + " end\n", + "\n", + " if size(J,1) == 1\n", + " W = one(J) - h*d*J\n", + " else\n", + " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", + " # M * dy/dt = F(t,y)\n", + " # we can simply replace eye(J) by M in the following expression\n", + " # (see Sec. 5 in [SR97])\n", + "\n", + " W = lufact( eye(J) - h*d*J )\n", + " end\n", + "\n", + " # approximate time-derivative of F\n", + " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", + "\n", + " # modified Rosenbrock formula\n", + " k1 = W\\(F0 + T)\n", + " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", + " k2 = W\\(F1 - k1) + k1\n", + " ynew = y + h*k2\n", + " F2 = F(t + h, ynew)\n", + " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", + "\n", + " err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate\n", + " delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error\n", + "\n", + " # check if new solution is acceptable\n", + " if err <= delta\n", + "\n", + " if points==:specified || points==:all\n", + " # only points in tspan are requested\n", + " # -> find relevant points in (t,t+h]\n", + " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", + " # rescale to (0,1]\n", + " s = (toi-t)/h\n", + "\n", + " # use interpolation formula to get solutions at t=toi\n", + " push!(tout, toi)\n", + " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", + " end\n", + " end\n", + " if (points==:all) && (tout[end]!=t+h)\n", + " # add the intermediate points\n", + " push!(tout, t + h)\n", + " push!(yout, ynew)\n", + " end\n", + "\n", + " # update solution\n", + " t = t + h\n", + " y = ynew\n", + "\n", + " F0 = F2 # use FSAL property\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + " # for new solution\n", + " end\n", + "\n", + " # update of the step size\n", + " h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) )\n", + " end\n", + "\n", + " return tout, yout\n", + "end\n", + "\n", + "\n", + "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", + "# with provided coefficients.\n", + "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", + "\n", + " if typeof(jacobian) == Function\n", + " G = jacobian\n", + " else\n", + " G = (t, x)->fdjacobian(F, x, t)\n", + " end\n", + "\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " solstep = 1\n", + " while solstep < length(tspan)\n", + " ts = tspan[solstep]\n", + " hs = h[solstep]\n", + " xs = x[solstep]\n", + " dFdx = G(ts, xs)\n", + " # FIXME\n", + " if size(dFdx,1) == 1\n", + " jac = 1/gamma/hs - dFdx[1]\n", + " else\n", + " jac = eye(dFdx)/gamma/hs - dFdx\n", + " end\n", + "\n", + " g = Array(typeof(x0), size(a,1))\n", + " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", + " x[solstep+1] = x[solstep] + b[1]*g[1]\n", + "\n", + " for i = 2:size(a,1)\n", + " dx = zero(x0)\n", + " dF = zero(x0/hs)\n", + " for j = 1:i-1\n", + " dx += a[i,j]*g[j]\n", + " dF += c[i,j]*g[j]\n", + " end\n", + " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", + " x[solstep+1] += b[i]*g[i]\n", + " end\n", + " solstep += 1\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "\n", + "# Kaps-Rentrop coefficients\n", + "const kr4_coefficients = (0.231,\n", + " [0 0 0 0\n", + " 2 0 0 0\n", + " 4.452470820736 4.16352878860 0 0\n", + " 4.452470820736 4.16352878860 0 0],\n", + " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", + " [ 0 0 0 0\n", + " -5.07167533877 0 0 0\n", + " 6.02015272865 0.1597500684673 0 0\n", + " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", + "\n", + "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Shampine coefficients\n", + "const s4_coefficients = (0.5,\n", + " [ 0 0 0 0\n", + " 2 0 0 0\n", + " 48/25 6/25 0 0\n", + " 48/25 6/25 0 0],\n", + " [19/9 1/2 25/108 125/108],\n", + " [ 0 0 0 0\n", + " -8 0 0 0\n", + " 372/25 12/5 0 0\n", + " -112/125 -54/125 -2/5 0],)\n", + "\n", + "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Use Shampine coefficients by default (matching Numerical Recipes)\n", + "const ode4s = ode4s_s\n", + "\n", + "const ms_coefficients4 = [ 1 0 0 0\n", + " -1/2 3/2 0 0\n", + " 5/12 -4/3 23/12 0\n", + " -9/24 37/24 -59/24 55/24]\n", + "\n", + "\n", + "#end # module ODE\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.2", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb new file mode 100644 index 000000000..8da022fb6 --- /dev/null +++ b/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb @@ -0,0 +1,582 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#isdefined(Base, :__precompile__) && __precompile__()\n", + "# Ordinary Differential Equation Solvers\n", + "\n", + "#module ODE\n", + "\n", + "using Polynomials\n", + "using Compat\n", + "\n", + "## minimal function export list\n", + "# adaptive non-stiff:\n", + "export ode23, ode45, ode78\n", + "# non-adaptive non-stiff:\n", + "export ode4, ode4ms\n", + "# adaptive stiff:\n", + "export ode23s\n", + "# non-adaptive stiff:\n", + "export ode4s\n", + "\n", + "## complete function export list: see runtests.jl\n", + "\n", + "###############################################################################\n", + "## Coefficient Tableaus\n", + "###############################################################################\n", + "\n", + "# Butcher Tableaus, or more generally coefficient tables\n", + "# see Hairer & Wanner 1992, p. 134, 166\n", + "\n", + "abstract Tableau{Name, S, T<:Real}\n", + "# Name is the name of the tableau/method (a symbol)\n", + "# S is the number of stages (an int)\n", + "# T is the type of the coefficients\n", + "#\n", + "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", + "#\n", + "# For all types of tableaus it assumes fields:\n", + "# order::(Int...) # order of the method(s)\n", + "#\n", + "# For Runge-Kutta methods it assumes fields:\n", + "# a::Matrix{T} # SxS matrix\n", + "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", + "# c::Vector{T} # S\n", + "#\n", + "# For a tableau:\n", + "# c1 | a_11 .... a_1s\n", + "# . | a_21 . .\n", + "# . | a_31 . .\n", + "# . | .... . .\n", + "# c_s | a_s1 ....... a_ss\n", + "# -----+--------------------\n", + "# | b_1 ... b_s this is the one used for stepping\n", + "# | b'_1 ... b'_s this is the one used for error-checking\n", + "\n", + "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", + "order(b::Tableau) = b.order\n", + "# Subtypes need to define a convert method to convert to a different\n", + "# eltype with signature:\n", + "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", + "\n", + "###############################################################################\n", + "## HELPER FUNCTIONS\n", + "###############################################################################\n", + "\n", + "# estimator for initial step based on book\n", + "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", + "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", + " # Returns first step, direction of integration and F evaluated at t0\n", + " tdir = sign(tend-t0)\n", + " tdir==0 && error(\"Zero time span\")\n", + " tau = max(reltol.*abs(x0), abstol)\n", + " d0 = norm(x0./tau, Inf)\n", + " f0 = F(t0, x0)\n", + " d1 = norm(f0./tau, Inf)\n", + " if d0 < 1e-5 || d1 < 1e-5\n", + " h0 = 1e-6\n", + " else\n", + " h0 = 0.01*(d0/d1)\n", + " end\n", + " # perform Euler step\n", + " x1 = x0 + tdir*h0*f0\n", + " f1 = F(t0 + tdir*h0, x1)\n", + " # estimate second derivative\n", + " d2 = norm((f1 - f0)./tau, Inf)/h0\n", + " if max(d1, d2) <= 1e-15\n", + " h1 = max(1e-6, 1e-3*h0)\n", + " else\n", + " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", + " h1 = 10.^pow\n", + " end\n", + " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", + "end\n", + "\n", + "# isoutofdomain takes the state and returns true if state is outside\n", + "# of the allowed domain. Used in adaptive step-control.\n", + "isoutofdomain(x) = isnan(x)\n", + "\n", + "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", + " # There are a few types involved in a call to a ODE solver which\n", + " # somehow need to be consistent:\n", + " #\n", + " # Et = eltype(tspan)\n", + " # Ey = eltype(y0)\n", + " # Ef = eltype(Tf)\n", + " #\n", + " # There are also the types of the containers, but they are not\n", + " # needed as `similar` is used to make containers.\n", + " # Tt = typeof(tspan)\n", + " # Ty = typeof(y0) # note, this can be a scalar\n", + " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", + " #\n", + " # Returns\n", + " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", + " # the moment a AbstractFloat\n", + " # - Eyf: suitable eltype of y and f(t,y)\n", + " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", + " # - Ty: container type of y0\n", + " # - btab: tableau with entries converted to Et\n", + "\n", + " # Needed interface:\n", + " # On components: /, -\n", + " # On container: eltype, promote_type\n", + " # On time container: eltype\n", + "\n", + " Ty = typeof(y0)\n", + " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", + "\n", + " Et = eltype(tspan)\n", + " @assert Et<:Real\n", + " if !(Et<:AbstractFloat)\n", + " Et = promote_type(Et, Float64)\n", + " end\n", + "\n", + " # if all are Floats, make them the same\n", + " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", + " Et = promote_type(Et, Eyf)\n", + " Eyf = Et\n", + " end\n", + "\n", + " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", + " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", + "\n", + " btab_ = convert(Et, btab)\n", + " return Et, Eyf, Ty, btab_\n", + "end\n", + "\n", + "###############################################################################\n", + "## NON-STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "include(\"runge_kutta.jl\")\n", + "\n", + "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", + "# with Adams-Bashforth-Moulton coefficients\n", + "function ode_ms(F, x0, tspan, order::Integer)\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " if 1 <= order <= 4\n", + " b = ms_coefficients4\n", + " else\n", + " b = zeros(order, order)\n", + " b[1:4, 1:4] = ms_coefficients4\n", + " for s = 5:order\n", + " for j = 0:(s - 1)\n", + " # Assign in correct order for multiplication below\n", + " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", + " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", + " b[s, s - j] = ((-1)^j / factorial(j)\n", + " / factorial(s - 1 - j) * polyval(p_int, 1))\n", + " end\n", + " end\n", + " end\n", + "\n", + " # TODO: use a better data structure here (should be an order-element circ buffer)\n", + " xdot = similar(x)\n", + " for i = 1:length(tspan)-1\n", + " # Need to run the first several steps at reduced order\n", + " steporder = min(i, order)\n", + " xdot[i] = F(tspan[i], x[i])\n", + "\n", + " x[i+1] = x[i]\n", + " for j = 1:steporder\n", + " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", + " end\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "# Use order 4 by default\n", + "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", + "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", + "\n", + "###############################################################################\n", + "## STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "# Crude forward finite differences estimator of Jacobian as fallback\n", + "\n", + "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", + "function fdjacobian(F, x::Number, t)\n", + " ftx = F(t, x)\n", + "\n", + " # The 100 below is heuristic\n", + " dx = (x .+ (x==0))./100\n", + " dFdx = (F(t,x+dx)-ftx)./dx\n", + "\n", + " return dFdx\n", + "end\n", + "\n", + "function fdjacobian(F, x, t)\n", + " ftx = F(t, x)\n", + " lx = max(length(x),1)\n", + " dFdx = zeros(eltype(x), lx, lx)\n", + " for j = 1:lx\n", + " # The 100 below is heuristic\n", + " dx = zeros(eltype(x), lx)\n", + " dx[j] = (x[j] .+ (x[j]==0))./100\n", + " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", + " end\n", + " return dFdx\n", + "end\n", + "\n", + "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", + "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", + "#\n", + "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", + "#\n", + "# supports keywords: points = :all | :specified (using dense output)\n", + "# jacobian = G(t,y)::Function | nothing (FD)\n", + "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", + " jacobian=nothing,\n", + " points=:all,\n", + " norm=Base.norm,\n", + " minstep=abs(tspan[end] - tspan[1])/1e18,\n", + " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", + " initstep=0.)\n", + "\n", + "\n", + " # select method for computing the Jacobian\n", + " if typeof(jacobian) == Function\n", + " jac = jacobian\n", + " else\n", + " # fallback finite-difference\n", + " jac = (t, y)->fdjacobian(F, y, t)\n", + " end\n", + "\n", + " # constants\n", + " const d = 1/(2 + sqrt(2))\n", + " const e32 = 6 + sqrt(2)\n", + "\n", + " # Check if norm provided by the user\n", + " normControl = (norm == Base.norm)\n", + " \n", + " @assert length(reltol) == 1 \"Relative tolerance must be a scalar\"\n", + " @assert length(abstol) == 1 || length(abstol) == length(y0) \"Dimension of Absolute tolerance does not match the dimension of the problem\"\n", + " if !normControl\n", + " @assert length(abstol) == 1 \"Norm control supports only scalar Absolute tolerance\"\n", + " end\n", + " \n", + " # Broadcast the abstol to a vector\n", + " if length(abstol) == 1 && length(y0) != 1 && normControl\n", + " abstol = abstol*ones(y0);\n", + " end\n", + " \n", + " # initialization\n", + " t = tspan[1]\n", + "\n", + " tfinal = tspan[end]\n", + "\n", + " h = initstep\n", + " if h == 0.\n", + " # initial guess at a step size\n", + " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", + " else\n", + " tdir = sign(tfinal - t)\n", + " F0 = F(t,y0)\n", + " end\n", + " h = tdir * min(abs(h), maxstep)\n", + "\n", + " y = y0\n", + " tout = Array(typeof(t), 1)\n", + " tout[1] = t # first output time\n", + " yout = Array(typeof(y0), 1)\n", + " yout[1] = deepcopy(y) # first output solution\n", + "\n", + "\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + "\n", + " while abs(t - tfinal) > 0 && minstep < abs(h)\n", + " if abs(t-tfinal) < abs(h)\n", + " h = tfinal - t\n", + " end\n", + "\n", + " if size(J,1) == 1\n", + " W = one(J) - h*d*J\n", + " else\n", + " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", + " # M * dy/dt = F(t,y)\n", + " # we can simply replace eye(J) by M in the following expression\n", + " # (see Sec. 5 in [SR97])\n", + "\n", + " W = lufact( eye(J) - h*d*J )\n", + " end\n", + "\n", + " # approximate time-derivative of F\n", + " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", + "\n", + " # modified Rosenbrock formula\n", + " k1 = W\\(F0 + T)\n", + " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", + " k2 = W\\(F1 - k1) + k1\n", + " ynew = y + h*k2\n", + " F2 = F(t + h, ynew)\n", + " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", + "\n", + " \n", + " # Component-wise error check similar to MATLAB ode23s\n", + " \n", + " threshold = abstol/reltol\n", + " \n", + " if normControl\n", + " err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate\n", + " else\n", + " err = (abs(h)/6)*(norm(k1 - 2*k2 + k3)/max(max(norm(y),norm(ynew)),threshold)) # error estimate\n", + " end\n", + " \n", + " # check if new solution is acceptable\n", + " if err <= reltol\n", + "\n", + " if points==:specified || points==:all\n", + " # only points in tspan are requested\n", + " # -> find relevant points in (t,t+h]\n", + " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", + " # rescale to (0,1]\n", + " s = (toi-t)/h\n", + "\n", + " # use interpolation formula to get solutions at t=toi\n", + " push!(tout, toi)\n", + " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", + " end\n", + " end\n", + " if (points==:all) && (tout[end]!=t+h)\n", + " # add the intermediate points\n", + " push!(tout, t + h)\n", + " push!(yout, ynew)\n", + " end\n", + "\n", + " # update solution\n", + " t = t + h\n", + " y = ynew\n", + "\n", + " F0 = F2 # use FSAL property\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + " # for new solution\n", + " end\n", + "\n", + " # update of the step size\n", + " h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) )\n", + " end\n", + "\n", + " return tout, yout\n", + "end\n", + "\n", + "\n", + "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", + "# with provided coefficients.\n", + "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", + "\n", + " if typeof(jacobian) == Function\n", + " G = jacobian\n", + " else\n", + " G = (t, x)->fdjacobian(F, x, t)\n", + " end\n", + "\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " solstep = 1\n", + " while solstep < length(tspan)\n", + " ts = tspan[solstep]\n", + " hs = h[solstep]\n", + " xs = x[solstep]\n", + " dFdx = G(ts, xs)\n", + " # FIXME\n", + " if size(dFdx,1) == 1\n", + " jac = 1/gamma/hs - dFdx[1]\n", + " else\n", + " jac = eye(dFdx)/gamma/hs - dFdx\n", + " end\n", + "\n", + " g = Array(typeof(x0), size(a,1))\n", + " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", + " x[solstep+1] = x[solstep] + b[1]*g[1]\n", + "\n", + " for i = 2:size(a,1)\n", + " dx = zero(x0)\n", + " dF = zero(x0/hs)\n", + " for j = 1:i-1\n", + " dx += a[i,j]*g[j]\n", + " dF += c[i,j]*g[j]\n", + " end\n", + " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", + " x[solstep+1] += b[i]*g[i]\n", + " end\n", + " solstep += 1\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "\n", + "# Kaps-Rentrop coefficients\n", + "const kr4_coefficients = (0.231,\n", + " [0 0 0 0\n", + " 2 0 0 0\n", + " 4.452470820736 4.16352878860 0 0\n", + " 4.452470820736 4.16352878860 0 0],\n", + " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", + " [ 0 0 0 0\n", + " -5.07167533877 0 0 0\n", + " 6.02015272865 0.1597500684673 0 0\n", + " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", + "\n", + "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Shampine coefficients\n", + "const s4_coefficients = (0.5,\n", + " [ 0 0 0 0\n", + " 2 0 0 0\n", + " 48/25 6/25 0 0\n", + " 48/25 6/25 0 0],\n", + " [19/9 1/2 25/108 125/108],\n", + " [ 0 0 0 0\n", + " -8 0 0 0\n", + " 372/25 12/5 0 0\n", + " -112/125 -54/125 -2/5 0],)\n", + "\n", + "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Use Shampine coefficients by default (matching Numerical Recipes)\n", + "const ode4s = ode4s_s\n", + "\n", + "const ms_coefficients4 = [ 1 0 0 0\n", + " -1/2 3/2 0 0\n", + " 5/12 -4/3 23/12 0\n", + " -9/24 37/24 -59/24 55/24]\n", + "\n", + "\n", + "#end # module ODE\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "tol = 1e-2\n", + "function norm1(x,p=1)\n", + " return norm(x,p)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;])\n", + "ys = hcat(y...).'\n", + "maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->y, 1., [1:-.001:0;])\n", + "maximum(abs(y-e.^(t-1))) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->y, 1., [0:.001:1;])\n", + "maximum(abs(y-e.^t)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->2t, 0., [0:.001:1;])\n", + "maximum(abs(y-t.^2)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->6.0, 0., [0:.1:1;];norm=norm1)\n", + "maximum(abs(y-6t)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "methods(norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.2", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/ODE.jl b/src/ODE.jl index e8ccc55c1..52374c8e9 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -66,10 +66,10 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) # Returns first step, direction of integration and F evaluated at t0 tdir = sign(tend-t0) tdir==0 && error("Zero time span") - tau = max(reltol*norm(x0, Inf), abstol) - d0 = norm(x0, Inf)/tau + tau = max(reltol.*abs(x0), abstol) + d0 = norm(x0./tau, Inf) f0 = F(t0, x0) - d1 = norm(f0, Inf)/tau + d1 = norm(f0./tau, Inf) if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 else @@ -79,7 +79,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) x1 = x0 + tdir*h0*f0 f1 = F(t0 + tdir*h0, x1) # estimate second derivative - d2 = norm(f1 - f0, Inf)/(tau*h0) + d2 = norm((f1 - f0)./tau, Inf)/h0 if max(d1, d2) <= 1e-15 h1 = max(1e-6, 1e-3*h0) else @@ -247,8 +247,15 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # constants const d = 1/(2 + sqrt(2)) const e32 = 6 + sqrt(2) - - + + @assert length(reltol) == 1 "Relative tolerance must be a scalar" + @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" + + # Broadcast the abstol to a vector + if length(abstol) == 1 && length(y0) != 1 + abstol = abstol*ones(y0); + end + # initialization t = tspan[1] @@ -300,11 +307,14 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, F2 = F(t + h, ynew) k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) - err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate - delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error - + + # Component-wise error check similar to MATLAB ode23s + + threshold = abstol/reltol + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate + # check if new solution is acceptable - if err <= delta + if err <= reltol if points==:specified || points==:all # only points in tspan are requested @@ -334,7 +344,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, end # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) ) + h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) ) end return tout, yout diff --git a/src/Untitled.ipynb b/src/Untitled.ipynb new file mode 100644 index 000000000..8da022fb6 --- /dev/null +++ b/src/Untitled.ipynb @@ -0,0 +1,582 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#isdefined(Base, :__precompile__) && __precompile__()\n", + "# Ordinary Differential Equation Solvers\n", + "\n", + "#module ODE\n", + "\n", + "using Polynomials\n", + "using Compat\n", + "\n", + "## minimal function export list\n", + "# adaptive non-stiff:\n", + "export ode23, ode45, ode78\n", + "# non-adaptive non-stiff:\n", + "export ode4, ode4ms\n", + "# adaptive stiff:\n", + "export ode23s\n", + "# non-adaptive stiff:\n", + "export ode4s\n", + "\n", + "## complete function export list: see runtests.jl\n", + "\n", + "###############################################################################\n", + "## Coefficient Tableaus\n", + "###############################################################################\n", + "\n", + "# Butcher Tableaus, or more generally coefficient tables\n", + "# see Hairer & Wanner 1992, p. 134, 166\n", + "\n", + "abstract Tableau{Name, S, T<:Real}\n", + "# Name is the name of the tableau/method (a symbol)\n", + "# S is the number of stages (an int)\n", + "# T is the type of the coefficients\n", + "#\n", + "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", + "#\n", + "# For all types of tableaus it assumes fields:\n", + "# order::(Int...) # order of the method(s)\n", + "#\n", + "# For Runge-Kutta methods it assumes fields:\n", + "# a::Matrix{T} # SxS matrix\n", + "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", + "# c::Vector{T} # S\n", + "#\n", + "# For a tableau:\n", + "# c1 | a_11 .... a_1s\n", + "# . | a_21 . .\n", + "# . | a_31 . .\n", + "# . | .... . .\n", + "# c_s | a_s1 ....... a_ss\n", + "# -----+--------------------\n", + "# | b_1 ... b_s this is the one used for stepping\n", + "# | b'_1 ... b'_s this is the one used for error-checking\n", + "\n", + "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", + "order(b::Tableau) = b.order\n", + "# Subtypes need to define a convert method to convert to a different\n", + "# eltype with signature:\n", + "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", + "\n", + "###############################################################################\n", + "## HELPER FUNCTIONS\n", + "###############################################################################\n", + "\n", + "# estimator for initial step based on book\n", + "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", + "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", + " # Returns first step, direction of integration and F evaluated at t0\n", + " tdir = sign(tend-t0)\n", + " tdir==0 && error(\"Zero time span\")\n", + " tau = max(reltol.*abs(x0), abstol)\n", + " d0 = norm(x0./tau, Inf)\n", + " f0 = F(t0, x0)\n", + " d1 = norm(f0./tau, Inf)\n", + " if d0 < 1e-5 || d1 < 1e-5\n", + " h0 = 1e-6\n", + " else\n", + " h0 = 0.01*(d0/d1)\n", + " end\n", + " # perform Euler step\n", + " x1 = x0 + tdir*h0*f0\n", + " f1 = F(t0 + tdir*h0, x1)\n", + " # estimate second derivative\n", + " d2 = norm((f1 - f0)./tau, Inf)/h0\n", + " if max(d1, d2) <= 1e-15\n", + " h1 = max(1e-6, 1e-3*h0)\n", + " else\n", + " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", + " h1 = 10.^pow\n", + " end\n", + " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", + "end\n", + "\n", + "# isoutofdomain takes the state and returns true if state is outside\n", + "# of the allowed domain. Used in adaptive step-control.\n", + "isoutofdomain(x) = isnan(x)\n", + "\n", + "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", + " # There are a few types involved in a call to a ODE solver which\n", + " # somehow need to be consistent:\n", + " #\n", + " # Et = eltype(tspan)\n", + " # Ey = eltype(y0)\n", + " # Ef = eltype(Tf)\n", + " #\n", + " # There are also the types of the containers, but they are not\n", + " # needed as `similar` is used to make containers.\n", + " # Tt = typeof(tspan)\n", + " # Ty = typeof(y0) # note, this can be a scalar\n", + " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", + " #\n", + " # Returns\n", + " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", + " # the moment a AbstractFloat\n", + " # - Eyf: suitable eltype of y and f(t,y)\n", + " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", + " # - Ty: container type of y0\n", + " # - btab: tableau with entries converted to Et\n", + "\n", + " # Needed interface:\n", + " # On components: /, -\n", + " # On container: eltype, promote_type\n", + " # On time container: eltype\n", + "\n", + " Ty = typeof(y0)\n", + " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", + "\n", + " Et = eltype(tspan)\n", + " @assert Et<:Real\n", + " if !(Et<:AbstractFloat)\n", + " Et = promote_type(Et, Float64)\n", + " end\n", + "\n", + " # if all are Floats, make them the same\n", + " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", + " Et = promote_type(Et, Eyf)\n", + " Eyf = Et\n", + " end\n", + "\n", + " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", + " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", + "\n", + " btab_ = convert(Et, btab)\n", + " return Et, Eyf, Ty, btab_\n", + "end\n", + "\n", + "###############################################################################\n", + "## NON-STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "include(\"runge_kutta.jl\")\n", + "\n", + "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", + "# with Adams-Bashforth-Moulton coefficients\n", + "function ode_ms(F, x0, tspan, order::Integer)\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " if 1 <= order <= 4\n", + " b = ms_coefficients4\n", + " else\n", + " b = zeros(order, order)\n", + " b[1:4, 1:4] = ms_coefficients4\n", + " for s = 5:order\n", + " for j = 0:(s - 1)\n", + " # Assign in correct order for multiplication below\n", + " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", + " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", + " b[s, s - j] = ((-1)^j / factorial(j)\n", + " / factorial(s - 1 - j) * polyval(p_int, 1))\n", + " end\n", + " end\n", + " end\n", + "\n", + " # TODO: use a better data structure here (should be an order-element circ buffer)\n", + " xdot = similar(x)\n", + " for i = 1:length(tspan)-1\n", + " # Need to run the first several steps at reduced order\n", + " steporder = min(i, order)\n", + " xdot[i] = F(tspan[i], x[i])\n", + "\n", + " x[i+1] = x[i]\n", + " for j = 1:steporder\n", + " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", + " end\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "# Use order 4 by default\n", + "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", + "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", + "\n", + "###############################################################################\n", + "## STIFF SOLVERS\n", + "###############################################################################\n", + "\n", + "# Crude forward finite differences estimator of Jacobian as fallback\n", + "\n", + "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", + "function fdjacobian(F, x::Number, t)\n", + " ftx = F(t, x)\n", + "\n", + " # The 100 below is heuristic\n", + " dx = (x .+ (x==0))./100\n", + " dFdx = (F(t,x+dx)-ftx)./dx\n", + "\n", + " return dFdx\n", + "end\n", + "\n", + "function fdjacobian(F, x, t)\n", + " ftx = F(t, x)\n", + " lx = max(length(x),1)\n", + " dFdx = zeros(eltype(x), lx, lx)\n", + " for j = 1:lx\n", + " # The 100 below is heuristic\n", + " dx = zeros(eltype(x), lx)\n", + " dx[j] = (x[j] .+ (x[j]==0))./100\n", + " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", + " end\n", + " return dFdx\n", + "end\n", + "\n", + "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", + "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", + "#\n", + "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", + "#\n", + "# supports keywords: points = :all | :specified (using dense output)\n", + "# jacobian = G(t,y)::Function | nothing (FD)\n", + "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", + " jacobian=nothing,\n", + " points=:all,\n", + " norm=Base.norm,\n", + " minstep=abs(tspan[end] - tspan[1])/1e18,\n", + " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", + " initstep=0.)\n", + "\n", + "\n", + " # select method for computing the Jacobian\n", + " if typeof(jacobian) == Function\n", + " jac = jacobian\n", + " else\n", + " # fallback finite-difference\n", + " jac = (t, y)->fdjacobian(F, y, t)\n", + " end\n", + "\n", + " # constants\n", + " const d = 1/(2 + sqrt(2))\n", + " const e32 = 6 + sqrt(2)\n", + "\n", + " # Check if norm provided by the user\n", + " normControl = (norm == Base.norm)\n", + " \n", + " @assert length(reltol) == 1 \"Relative tolerance must be a scalar\"\n", + " @assert length(abstol) == 1 || length(abstol) == length(y0) \"Dimension of Absolute tolerance does not match the dimension of the problem\"\n", + " if !normControl\n", + " @assert length(abstol) == 1 \"Norm control supports only scalar Absolute tolerance\"\n", + " end\n", + " \n", + " # Broadcast the abstol to a vector\n", + " if length(abstol) == 1 && length(y0) != 1 && normControl\n", + " abstol = abstol*ones(y0);\n", + " end\n", + " \n", + " # initialization\n", + " t = tspan[1]\n", + "\n", + " tfinal = tspan[end]\n", + "\n", + " h = initstep\n", + " if h == 0.\n", + " # initial guess at a step size\n", + " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", + " else\n", + " tdir = sign(tfinal - t)\n", + " F0 = F(t,y0)\n", + " end\n", + " h = tdir * min(abs(h), maxstep)\n", + "\n", + " y = y0\n", + " tout = Array(typeof(t), 1)\n", + " tout[1] = t # first output time\n", + " yout = Array(typeof(y0), 1)\n", + " yout[1] = deepcopy(y) # first output solution\n", + "\n", + "\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + "\n", + " while abs(t - tfinal) > 0 && minstep < abs(h)\n", + " if abs(t-tfinal) < abs(h)\n", + " h = tfinal - t\n", + " end\n", + "\n", + " if size(J,1) == 1\n", + " W = one(J) - h*d*J\n", + " else\n", + " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", + " # M * dy/dt = F(t,y)\n", + " # we can simply replace eye(J) by M in the following expression\n", + " # (see Sec. 5 in [SR97])\n", + "\n", + " W = lufact( eye(J) - h*d*J )\n", + " end\n", + "\n", + " # approximate time-derivative of F\n", + " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", + "\n", + " # modified Rosenbrock formula\n", + " k1 = W\\(F0 + T)\n", + " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", + " k2 = W\\(F1 - k1) + k1\n", + " ynew = y + h*k2\n", + " F2 = F(t + h, ynew)\n", + " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", + "\n", + " \n", + " # Component-wise error check similar to MATLAB ode23s\n", + " \n", + " threshold = abstol/reltol\n", + " \n", + " if normControl\n", + " err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate\n", + " else\n", + " err = (abs(h)/6)*(norm(k1 - 2*k2 + k3)/max(max(norm(y),norm(ynew)),threshold)) # error estimate\n", + " end\n", + " \n", + " # check if new solution is acceptable\n", + " if err <= reltol\n", + "\n", + " if points==:specified || points==:all\n", + " # only points in tspan are requested\n", + " # -> find relevant points in (t,t+h]\n", + " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", + " # rescale to (0,1]\n", + " s = (toi-t)/h\n", + "\n", + " # use interpolation formula to get solutions at t=toi\n", + " push!(tout, toi)\n", + " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", + " end\n", + " end\n", + " if (points==:all) && (tout[end]!=t+h)\n", + " # add the intermediate points\n", + " push!(tout, t + h)\n", + " push!(yout, ynew)\n", + " end\n", + "\n", + " # update solution\n", + " t = t + h\n", + " y = ynew\n", + "\n", + " F0 = F2 # use FSAL property\n", + " J = jac(t,y) # get Jacobian of F wrt y\n", + " # for new solution\n", + " end\n", + "\n", + " # update of the step size\n", + " h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) )\n", + " end\n", + "\n", + " return tout, yout\n", + "end\n", + "\n", + "\n", + "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", + "# with provided coefficients.\n", + "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", + "\n", + " if typeof(jacobian) == Function\n", + " G = jacobian\n", + " else\n", + " G = (t, x)->fdjacobian(F, x, t)\n", + " end\n", + "\n", + " h = diff(tspan)\n", + " x = Array(typeof(x0), length(tspan))\n", + " x[1] = x0\n", + "\n", + " solstep = 1\n", + " while solstep < length(tspan)\n", + " ts = tspan[solstep]\n", + " hs = h[solstep]\n", + " xs = x[solstep]\n", + " dFdx = G(ts, xs)\n", + " # FIXME\n", + " if size(dFdx,1) == 1\n", + " jac = 1/gamma/hs - dFdx[1]\n", + " else\n", + " jac = eye(dFdx)/gamma/hs - dFdx\n", + " end\n", + "\n", + " g = Array(typeof(x0), size(a,1))\n", + " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", + " x[solstep+1] = x[solstep] + b[1]*g[1]\n", + "\n", + " for i = 2:size(a,1)\n", + " dx = zero(x0)\n", + " dF = zero(x0/hs)\n", + " for j = 1:i-1\n", + " dx += a[i,j]*g[j]\n", + " dF += c[i,j]*g[j]\n", + " end\n", + " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", + " x[solstep+1] += b[i]*g[i]\n", + " end\n", + " solstep += 1\n", + " end\n", + " return vcat(tspan), x\n", + "end\n", + "\n", + "\n", + "# Kaps-Rentrop coefficients\n", + "const kr4_coefficients = (0.231,\n", + " [0 0 0 0\n", + " 2 0 0 0\n", + " 4.452470820736 4.16352878860 0 0\n", + " 4.452470820736 4.16352878860 0 0],\n", + " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", + " [ 0 0 0 0\n", + " -5.07167533877 0 0 0\n", + " 6.02015272865 0.1597500684673 0 0\n", + " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", + "\n", + "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Shampine coefficients\n", + "const s4_coefficients = (0.5,\n", + " [ 0 0 0 0\n", + " 2 0 0 0\n", + " 48/25 6/25 0 0\n", + " 48/25 6/25 0 0],\n", + " [19/9 1/2 25/108 125/108],\n", + " [ 0 0 0 0\n", + " -8 0 0 0\n", + " 372/25 12/5 0 0\n", + " -112/125 -54/125 -2/5 0],)\n", + "\n", + "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", + "\n", + "# Use Shampine coefficients by default (matching Numerical Recipes)\n", + "const ode4s = ode4s_s\n", + "\n", + "const ms_coefficients4 = [ 1 0 0 0\n", + " -1/2 3/2 0 0\n", + " 5/12 -4/3 23/12 0\n", + " -9/24 37/24 -59/24 55/24]\n", + "\n", + "\n", + "#end # module ODE\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "tol = 1e-2\n", + "function norm1(x,p=1)\n", + " return norm(x,p)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;])\n", + "ys = hcat(y...).'\n", + "maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->y, 1., [1:-.001:0;])\n", + "maximum(abs(y-e.^(t-1))) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->y, 1., [0:.001:1;])\n", + "maximum(abs(y-e.^t)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->2t, 0., [0:.001:1;])\n", + "maximum(abs(y-t.^2)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t,y=ode23s((t,y)->6.0, 0., [0:.1:1;];norm=norm1)\n", + "maximum(abs(y-6t)) < tol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "methods(norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.2", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index a080602de..5f1a994a6 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -250,6 +250,13 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici t = tspan[1] tstart = tspan[1] tend = tspan[end] + + if length(reltol) == 1 + reltol = reltol*ones(y0); + end + if length(abstol) == 1 + abstol = abstol*ones(y0); + end # work arrays: y = similar(y0, Eyf, dof) # y at time t @@ -412,7 +419,7 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, for d=1:dof # if outside of domain (usually NaN) then make step size smaller by maximum isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol + max(norm(x0[d]), norm(xtrial[d]))*reltol) # Eq 4.10 + xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 end err = norm(xerr, 2) # Eq. 4.11 newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified diff --git a/test/runtests.jl b/test/runtests.jl index 9c88cd281..0be0559f2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -81,7 +81,7 @@ let refsol = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050] # reference solution at tspan[2] - @test norm(refsol-y[end], Inf) < 2e-10 + @test norm(refsol-y[end], Inf) < 3e-10 end include("interface-tests.jl") From 42f6b4550189de9def6993e770f49532ad5f9781 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 17:07:19 +0200 Subject: [PATCH 02/14] Pre .gitignore changes --- .ipynb_checkpoints/Untitled-checkpoint.ipynb | 487 ------------------- test/runtests.jl | 6 +- 2 files changed, 5 insertions(+), 488 deletions(-) delete mode 100644 .ipynb_checkpoints/Untitled-checkpoint.ipynb diff --git a/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/.ipynb_checkpoints/Untitled-checkpoint.ipynb deleted file mode 100644 index 97d2e1d78..000000000 --- a/.ipynb_checkpoints/Untitled-checkpoint.ipynb +++ /dev/null @@ -1,487 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: could not open file /mnt/526E44596E4437CD/GitRepo/ODE.jl/runge_kutta.jl\nwhile loading In[1], in expression starting on line 149", - "output_type": "error", - "traceback": [ - "LoadError: could not open file /mnt/526E44596E4437CD/GitRepo/ODE.jl/runge_kutta.jl\nwhile loading In[1], in expression starting on line 149", - "", - " in include at ./boot.jl:261", - " in include_from_node1 at ./loading.jl:304" - ] - } - ], - "source": [ - "#isdefined(Base, :__precompile__) && __precompile__()\n", - "# Ordinary Differential Equation Solvers\n", - "\n", - "#module ODE\n", - "\n", - "using Polynomials\n", - "using Compat\n", - "\n", - "## minimal function export list\n", - "# adaptive non-stiff:\n", - "export ode23, ode45, ode78\n", - "# non-adaptive non-stiff:\n", - "export ode4, ode4ms\n", - "# adaptive stiff:\n", - "export ode23s\n", - "# non-adaptive stiff:\n", - "export ode4s\n", - "\n", - "## complete function export list: see runtests.jl\n", - "\n", - "###############################################################################\n", - "## Coefficient Tableaus\n", - "###############################################################################\n", - "\n", - "# Butcher Tableaus, or more generally coefficient tables\n", - "# see Hairer & Wanner 1992, p. 134, 166\n", - "\n", - "abstract Tableau{Name, S, T<:Real}\n", - "# Name is the name of the tableau/method (a symbol)\n", - "# S is the number of stages (an int)\n", - "# T is the type of the coefficients\n", - "#\n", - "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", - "#\n", - "# For all types of tableaus it assumes fields:\n", - "# order::(Int...) # order of the method(s)\n", - "#\n", - "# For Runge-Kutta methods it assumes fields:\n", - "# a::Matrix{T} # SxS matrix\n", - "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", - "# c::Vector{T} # S\n", - "#\n", - "# For a tableau:\n", - "# c1 | a_11 .... a_1s\n", - "# . | a_21 . .\n", - "# . | a_31 . .\n", - "# . | .... . .\n", - "# c_s | a_s1 ....... a_ss\n", - "# -----+--------------------\n", - "# | b_1 ... b_s this is the one used for stepping\n", - "# | b'_1 ... b'_s this is the one used for error-checking\n", - "\n", - "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", - "order(b::Tableau) = b.order\n", - "# Subtypes need to define a convert method to convert to a different\n", - "# eltype with signature:\n", - "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", - "\n", - "###############################################################################\n", - "## HELPER FUNCTIONS\n", - "###############################################################################\n", - "\n", - "# estimator for initial step based on book\n", - "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", - "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", - " # Returns first step, direction of integration and F evaluated at t0\n", - " tdir = sign(tend-t0)\n", - " tdir==0 && error(\"Zero time span\")\n", - " tau = max(reltol*norm(x0, Inf), abstol)\n", - " d0 = norm(x0, Inf)/tau\n", - " f0 = F(t0, x0)\n", - " d1 = norm(f0, Inf)/tau\n", - " if d0 < 1e-5 || d1 < 1e-5\n", - " h0 = 1e-6\n", - " else\n", - " h0 = 0.01*(d0/d1)\n", - " end\n", - " # perform Euler step\n", - " x1 = x0 + tdir*h0*f0\n", - " f1 = F(t0 + tdir*h0, x1)\n", - " # estimate second derivative\n", - " d2 = norm(f1 - f0, Inf)/(tau*h0)\n", - " if max(d1, d2) <= 1e-15\n", - " h1 = max(1e-6, 1e-3*h0)\n", - " else\n", - " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", - " h1 = 10.^pow\n", - " end\n", - " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", - "end\n", - "\n", - "# isoutofdomain takes the state and returns true if state is outside\n", - "# of the allowed domain. Used in adaptive step-control.\n", - "isoutofdomain(x) = isnan(x)\n", - "\n", - "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", - " # There are a few types involved in a call to a ODE solver which\n", - " # somehow need to be consistent:\n", - " #\n", - " # Et = eltype(tspan)\n", - " # Ey = eltype(y0)\n", - " # Ef = eltype(Tf)\n", - " #\n", - " # There are also the types of the containers, but they are not\n", - " # needed as `similar` is used to make containers.\n", - " # Tt = typeof(tspan)\n", - " # Ty = typeof(y0) # note, this can be a scalar\n", - " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", - " #\n", - " # Returns\n", - " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", - " # the moment a AbstractFloat\n", - " # - Eyf: suitable eltype of y and f(t,y)\n", - " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", - " # - Ty: container type of y0\n", - " # - btab: tableau with entries converted to Et\n", - "\n", - " # Needed interface:\n", - " # On components: /, -\n", - " # On container: eltype, promote_type\n", - " # On time container: eltype\n", - "\n", - " Ty = typeof(y0)\n", - " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", - "\n", - " Et = eltype(tspan)\n", - " @assert Et<:Real\n", - " if !(Et<:AbstractFloat)\n", - " Et = promote_type(Et, Float64)\n", - " end\n", - "\n", - " # if all are Floats, make them the same\n", - " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", - " Et = promote_type(Et, Eyf)\n", - " Eyf = Et\n", - " end\n", - "\n", - " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", - " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", - "\n", - " btab_ = convert(Et, btab)\n", - " return Et, Eyf, Ty, btab_\n", - "end\n", - "\n", - "###############################################################################\n", - "## NON-STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "include(\"runge_kutta.jl\")\n", - "\n", - "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", - "# with Adams-Bashforth-Moulton coefficients\n", - "function ode_ms(F, x0, tspan, order::Integer)\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " if 1 <= order <= 4\n", - " b = ms_coefficients4\n", - " else\n", - " b = zeros(order, order)\n", - " b[1:4, 1:4] = ms_coefficients4\n", - " for s = 5:order\n", - " for j = 0:(s - 1)\n", - " # Assign in correct order for multiplication below\n", - " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", - " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", - " b[s, s - j] = ((-1)^j / factorial(j)\n", - " / factorial(s - 1 - j) * polyval(p_int, 1))\n", - " end\n", - " end\n", - " end\n", - "\n", - " # TODO: use a better data structure here (should be an order-element circ buffer)\n", - " xdot = similar(x)\n", - " for i = 1:length(tspan)-1\n", - " # Need to run the first several steps at reduced order\n", - " steporder = min(i, order)\n", - " xdot[i] = F(tspan[i], x[i])\n", - "\n", - " x[i+1] = x[i]\n", - " for j = 1:steporder\n", - " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", - " end\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "# Use order 4 by default\n", - "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", - "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", - "\n", - "###############################################################################\n", - "## STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "# Crude forward finite differences estimator of Jacobian as fallback\n", - "\n", - "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", - "function fdjacobian(F, x::Number, t)\n", - " ftx = F(t, x)\n", - "\n", - " # The 100 below is heuristic\n", - " dx = (x .+ (x==0))./100\n", - " dFdx = (F(t,x+dx)-ftx)./dx\n", - "\n", - " return dFdx\n", - "end\n", - "\n", - "function fdjacobian(F, x, t)\n", - " ftx = F(t, x)\n", - " lx = max(length(x),1)\n", - " dFdx = zeros(eltype(x), lx, lx)\n", - " for j = 1:lx\n", - " # The 100 below is heuristic\n", - " dx = zeros(eltype(x), lx)\n", - " dx[j] = (x[j] .+ (x[j]==0))./100\n", - " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", - " end\n", - " return dFdx\n", - "end\n", - "\n", - "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", - "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", - "#\n", - "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", - "#\n", - "# supports keywords: points = :all | :specified (using dense output)\n", - "# jacobian = G(t,y)::Function | nothing (FD)\n", - "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", - " jacobian=nothing,\n", - " points=:all,\n", - " norm=Base.norm,\n", - " minstep=abs(tspan[end] - tspan[1])/1e18,\n", - " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", - " initstep=0.)\n", - "\n", - "\n", - " # select method for computing the Jacobian\n", - " if typeof(jacobian) == Function\n", - " jac = jacobian\n", - " else\n", - " # fallback finite-difference\n", - " jac = (t, y)->fdjacobian(F, y, t)\n", - " end\n", - "\n", - " # constants\n", - " const d = 1/(2 + sqrt(2))\n", - " const e32 = 6 + sqrt(2)\n", - "\n", - "\n", - " # initialization\n", - " t = tspan[1]\n", - "\n", - " tfinal = tspan[end]\n", - "\n", - " h = initstep\n", - " if h == 0.\n", - " # initial guess at a step size\n", - " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", - " else\n", - " tdir = sign(tfinal - t)\n", - " F0 = F(t,y0)\n", - " end\n", - " h = tdir * min(abs(h), maxstep)\n", - "\n", - " y = y0\n", - " tout = Array(typeof(t), 1)\n", - " tout[1] = t # first output time\n", - " yout = Array(typeof(y0), 1)\n", - " yout[1] = deepcopy(y) # first output solution\n", - "\n", - "\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - "\n", - " while abs(t - tfinal) > 0 && minstep < abs(h)\n", - " if abs(t-tfinal) < abs(h)\n", - " h = tfinal - t\n", - " end\n", - "\n", - " if size(J,1) == 1\n", - " W = one(J) - h*d*J\n", - " else\n", - " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", - " # M * dy/dt = F(t,y)\n", - " # we can simply replace eye(J) by M in the following expression\n", - " # (see Sec. 5 in [SR97])\n", - "\n", - " W = lufact( eye(J) - h*d*J )\n", - " end\n", - "\n", - " # approximate time-derivative of F\n", - " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", - "\n", - " # modified Rosenbrock formula\n", - " k1 = W\\(F0 + T)\n", - " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", - " k2 = W\\(F1 - k1) + k1\n", - " ynew = y + h*k2\n", - " F2 = F(t + h, ynew)\n", - " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", - "\n", - " err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate\n", - " delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error\n", - "\n", - " # check if new solution is acceptable\n", - " if err <= delta\n", - "\n", - " if points==:specified || points==:all\n", - " # only points in tspan are requested\n", - " # -> find relevant points in (t,t+h]\n", - " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", - " # rescale to (0,1]\n", - " s = (toi-t)/h\n", - "\n", - " # use interpolation formula to get solutions at t=toi\n", - " push!(tout, toi)\n", - " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", - " end\n", - " end\n", - " if (points==:all) && (tout[end]!=t+h)\n", - " # add the intermediate points\n", - " push!(tout, t + h)\n", - " push!(yout, ynew)\n", - " end\n", - "\n", - " # update solution\n", - " t = t + h\n", - " y = ynew\n", - "\n", - " F0 = F2 # use FSAL property\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - " # for new solution\n", - " end\n", - "\n", - " # update of the step size\n", - " h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) )\n", - " end\n", - "\n", - " return tout, yout\n", - "end\n", - "\n", - "\n", - "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", - "# with provided coefficients.\n", - "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", - "\n", - " if typeof(jacobian) == Function\n", - " G = jacobian\n", - " else\n", - " G = (t, x)->fdjacobian(F, x, t)\n", - " end\n", - "\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " solstep = 1\n", - " while solstep < length(tspan)\n", - " ts = tspan[solstep]\n", - " hs = h[solstep]\n", - " xs = x[solstep]\n", - " dFdx = G(ts, xs)\n", - " # FIXME\n", - " if size(dFdx,1) == 1\n", - " jac = 1/gamma/hs - dFdx[1]\n", - " else\n", - " jac = eye(dFdx)/gamma/hs - dFdx\n", - " end\n", - "\n", - " g = Array(typeof(x0), size(a,1))\n", - " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", - " x[solstep+1] = x[solstep] + b[1]*g[1]\n", - "\n", - " for i = 2:size(a,1)\n", - " dx = zero(x0)\n", - " dF = zero(x0/hs)\n", - " for j = 1:i-1\n", - " dx += a[i,j]*g[j]\n", - " dF += c[i,j]*g[j]\n", - " end\n", - " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", - " x[solstep+1] += b[i]*g[i]\n", - " end\n", - " solstep += 1\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "\n", - "# Kaps-Rentrop coefficients\n", - "const kr4_coefficients = (0.231,\n", - " [0 0 0 0\n", - " 2 0 0 0\n", - " 4.452470820736 4.16352878860 0 0\n", - " 4.452470820736 4.16352878860 0 0],\n", - " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", - " [ 0 0 0 0\n", - " -5.07167533877 0 0 0\n", - " 6.02015272865 0.1597500684673 0 0\n", - " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", - "\n", - "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Shampine coefficients\n", - "const s4_coefficients = (0.5,\n", - " [ 0 0 0 0\n", - " 2 0 0 0\n", - " 48/25 6/25 0 0\n", - " 48/25 6/25 0 0],\n", - " [19/9 1/2 25/108 125/108],\n", - " [ 0 0 0 0\n", - " -8 0 0 0\n", - " 372/25 12/5 0 0\n", - " -112/125 -54/125 -2/5 0],)\n", - "\n", - "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Use Shampine coefficients by default (matching Numerical Recipes)\n", - "const ode4s = ode4s_s\n", - "\n", - "const ms_coefficients4 = [ 1 0 0 0\n", - " -1/2 3/2 0 0\n", - " 5/12 -4/3 23/12 0\n", - " -9/24 37/24 -59/24 55/24]\n", - "\n", - "\n", - "#end # module ODE\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 0.4.2", - "language": "julia", - "name": "julia-0.4" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "0.4.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} diff --git a/test/runtests.jl b/test/runtests.jl index 0be0559f2..c2c211178 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,10 @@ end # rober testcase from http://www.unige.ch/~hairer/testset/testset.html +# Changing the test value. +# Even the Matlab ode23s does not give the error to be less than 2e-10. +# The value comes out to be approximately 2.042354e-10. +# Hence changing the test value to 2.1e-10 let println("ROBER test case") function f(t, y) @@ -81,7 +85,7 @@ let refsol = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050] # reference solution at tspan[2] - @test norm(refsol-y[end], Inf) < 3e-10 + @test norm(refsol-y[end], Inf) < 2.1e-10 end include("interface-tests.jl") From 0eeb96bb532e4a97fd23a408029480cef7694865 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 17:08:22 +0200 Subject: [PATCH 03/14] Post .gitignore changes --- .travis.yml | 17 - .../Untitled-checkpoint.ipynb | 582 ------------------ 2 files changed, 599 deletions(-) delete mode 100644 .travis.yml delete mode 100644 src/.ipynb_checkpoints/Untitled-checkpoint.ipynb diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 1fc28491b..000000000 --- a/.travis.yml +++ /dev/null @@ -1,17 +0,0 @@ -language: julia -os: - - osx - - linux -julia: - - 0.3 - - release - - nightly -git: - depth: 999999 -notifications: - email: false -script: - - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi - - julia -e 'Pkg.clone(pwd()); Pkg.build("ODE"); Pkg.test("ODE"; coverage=true)'; -after_success: - - julia -e 'cd(Pkg.dir("ODE")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'; diff --git a/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb deleted file mode 100644 index 8da022fb6..000000000 --- a/src/.ipynb_checkpoints/Untitled-checkpoint.ipynb +++ /dev/null @@ -1,582 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "#isdefined(Base, :__precompile__) && __precompile__()\n", - "# Ordinary Differential Equation Solvers\n", - "\n", - "#module ODE\n", - "\n", - "using Polynomials\n", - "using Compat\n", - "\n", - "## minimal function export list\n", - "# adaptive non-stiff:\n", - "export ode23, ode45, ode78\n", - "# non-adaptive non-stiff:\n", - "export ode4, ode4ms\n", - "# adaptive stiff:\n", - "export ode23s\n", - "# non-adaptive stiff:\n", - "export ode4s\n", - "\n", - "## complete function export list: see runtests.jl\n", - "\n", - "###############################################################################\n", - "## Coefficient Tableaus\n", - "###############################################################################\n", - "\n", - "# Butcher Tableaus, or more generally coefficient tables\n", - "# see Hairer & Wanner 1992, p. 134, 166\n", - "\n", - "abstract Tableau{Name, S, T<:Real}\n", - "# Name is the name of the tableau/method (a symbol)\n", - "# S is the number of stages (an int)\n", - "# T is the type of the coefficients\n", - "#\n", - "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", - "#\n", - "# For all types of tableaus it assumes fields:\n", - "# order::(Int...) # order of the method(s)\n", - "#\n", - "# For Runge-Kutta methods it assumes fields:\n", - "# a::Matrix{T} # SxS matrix\n", - "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", - "# c::Vector{T} # S\n", - "#\n", - "# For a tableau:\n", - "# c1 | a_11 .... a_1s\n", - "# . | a_21 . .\n", - "# . | a_31 . .\n", - "# . | .... . .\n", - "# c_s | a_s1 ....... a_ss\n", - "# -----+--------------------\n", - "# | b_1 ... b_s this is the one used for stepping\n", - "# | b'_1 ... b'_s this is the one used for error-checking\n", - "\n", - "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", - "order(b::Tableau) = b.order\n", - "# Subtypes need to define a convert method to convert to a different\n", - "# eltype with signature:\n", - "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", - "\n", - "###############################################################################\n", - "## HELPER FUNCTIONS\n", - "###############################################################################\n", - "\n", - "# estimator for initial step based on book\n", - "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", - "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", - " # Returns first step, direction of integration and F evaluated at t0\n", - " tdir = sign(tend-t0)\n", - " tdir==0 && error(\"Zero time span\")\n", - " tau = max(reltol.*abs(x0), abstol)\n", - " d0 = norm(x0./tau, Inf)\n", - " f0 = F(t0, x0)\n", - " d1 = norm(f0./tau, Inf)\n", - " if d0 < 1e-5 || d1 < 1e-5\n", - " h0 = 1e-6\n", - " else\n", - " h0 = 0.01*(d0/d1)\n", - " end\n", - " # perform Euler step\n", - " x1 = x0 + tdir*h0*f0\n", - " f1 = F(t0 + tdir*h0, x1)\n", - " # estimate second derivative\n", - " d2 = norm((f1 - f0)./tau, Inf)/h0\n", - " if max(d1, d2) <= 1e-15\n", - " h1 = max(1e-6, 1e-3*h0)\n", - " else\n", - " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", - " h1 = 10.^pow\n", - " end\n", - " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", - "end\n", - "\n", - "# isoutofdomain takes the state and returns true if state is outside\n", - "# of the allowed domain. Used in adaptive step-control.\n", - "isoutofdomain(x) = isnan(x)\n", - "\n", - "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", - " # There are a few types involved in a call to a ODE solver which\n", - " # somehow need to be consistent:\n", - " #\n", - " # Et = eltype(tspan)\n", - " # Ey = eltype(y0)\n", - " # Ef = eltype(Tf)\n", - " #\n", - " # There are also the types of the containers, but they are not\n", - " # needed as `similar` is used to make containers.\n", - " # Tt = typeof(tspan)\n", - " # Ty = typeof(y0) # note, this can be a scalar\n", - " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", - " #\n", - " # Returns\n", - " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", - " # the moment a AbstractFloat\n", - " # - Eyf: suitable eltype of y and f(t,y)\n", - " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", - " # - Ty: container type of y0\n", - " # - btab: tableau with entries converted to Et\n", - "\n", - " # Needed interface:\n", - " # On components: /, -\n", - " # On container: eltype, promote_type\n", - " # On time container: eltype\n", - "\n", - " Ty = typeof(y0)\n", - " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", - "\n", - " Et = eltype(tspan)\n", - " @assert Et<:Real\n", - " if !(Et<:AbstractFloat)\n", - " Et = promote_type(Et, Float64)\n", - " end\n", - "\n", - " # if all are Floats, make them the same\n", - " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", - " Et = promote_type(Et, Eyf)\n", - " Eyf = Et\n", - " end\n", - "\n", - " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", - " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", - "\n", - " btab_ = convert(Et, btab)\n", - " return Et, Eyf, Ty, btab_\n", - "end\n", - "\n", - "###############################################################################\n", - "## NON-STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "include(\"runge_kutta.jl\")\n", - "\n", - "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", - "# with Adams-Bashforth-Moulton coefficients\n", - "function ode_ms(F, x0, tspan, order::Integer)\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " if 1 <= order <= 4\n", - " b = ms_coefficients4\n", - " else\n", - " b = zeros(order, order)\n", - " b[1:4, 1:4] = ms_coefficients4\n", - " for s = 5:order\n", - " for j = 0:(s - 1)\n", - " # Assign in correct order for multiplication below\n", - " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", - " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", - " b[s, s - j] = ((-1)^j / factorial(j)\n", - " / factorial(s - 1 - j) * polyval(p_int, 1))\n", - " end\n", - " end\n", - " end\n", - "\n", - " # TODO: use a better data structure here (should be an order-element circ buffer)\n", - " xdot = similar(x)\n", - " for i = 1:length(tspan)-1\n", - " # Need to run the first several steps at reduced order\n", - " steporder = min(i, order)\n", - " xdot[i] = F(tspan[i], x[i])\n", - "\n", - " x[i+1] = x[i]\n", - " for j = 1:steporder\n", - " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", - " end\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "# Use order 4 by default\n", - "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", - "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", - "\n", - "###############################################################################\n", - "## STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "# Crude forward finite differences estimator of Jacobian as fallback\n", - "\n", - "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", - "function fdjacobian(F, x::Number, t)\n", - " ftx = F(t, x)\n", - "\n", - " # The 100 below is heuristic\n", - " dx = (x .+ (x==0))./100\n", - " dFdx = (F(t,x+dx)-ftx)./dx\n", - "\n", - " return dFdx\n", - "end\n", - "\n", - "function fdjacobian(F, x, t)\n", - " ftx = F(t, x)\n", - " lx = max(length(x),1)\n", - " dFdx = zeros(eltype(x), lx, lx)\n", - " for j = 1:lx\n", - " # The 100 below is heuristic\n", - " dx = zeros(eltype(x), lx)\n", - " dx[j] = (x[j] .+ (x[j]==0))./100\n", - " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", - " end\n", - " return dFdx\n", - "end\n", - "\n", - "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", - "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", - "#\n", - "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", - "#\n", - "# supports keywords: points = :all | :specified (using dense output)\n", - "# jacobian = G(t,y)::Function | nothing (FD)\n", - "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", - " jacobian=nothing,\n", - " points=:all,\n", - " norm=Base.norm,\n", - " minstep=abs(tspan[end] - tspan[1])/1e18,\n", - " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", - " initstep=0.)\n", - "\n", - "\n", - " # select method for computing the Jacobian\n", - " if typeof(jacobian) == Function\n", - " jac = jacobian\n", - " else\n", - " # fallback finite-difference\n", - " jac = (t, y)->fdjacobian(F, y, t)\n", - " end\n", - "\n", - " # constants\n", - " const d = 1/(2 + sqrt(2))\n", - " const e32 = 6 + sqrt(2)\n", - "\n", - " # Check if norm provided by the user\n", - " normControl = (norm == Base.norm)\n", - " \n", - " @assert length(reltol) == 1 \"Relative tolerance must be a scalar\"\n", - " @assert length(abstol) == 1 || length(abstol) == length(y0) \"Dimension of Absolute tolerance does not match the dimension of the problem\"\n", - " if !normControl\n", - " @assert length(abstol) == 1 \"Norm control supports only scalar Absolute tolerance\"\n", - " end\n", - " \n", - " # Broadcast the abstol to a vector\n", - " if length(abstol) == 1 && length(y0) != 1 && normControl\n", - " abstol = abstol*ones(y0);\n", - " end\n", - " \n", - " # initialization\n", - " t = tspan[1]\n", - "\n", - " tfinal = tspan[end]\n", - "\n", - " h = initstep\n", - " if h == 0.\n", - " # initial guess at a step size\n", - " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", - " else\n", - " tdir = sign(tfinal - t)\n", - " F0 = F(t,y0)\n", - " end\n", - " h = tdir * min(abs(h), maxstep)\n", - "\n", - " y = y0\n", - " tout = Array(typeof(t), 1)\n", - " tout[1] = t # first output time\n", - " yout = Array(typeof(y0), 1)\n", - " yout[1] = deepcopy(y) # first output solution\n", - "\n", - "\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - "\n", - " while abs(t - tfinal) > 0 && minstep < abs(h)\n", - " if abs(t-tfinal) < abs(h)\n", - " h = tfinal - t\n", - " end\n", - "\n", - " if size(J,1) == 1\n", - " W = one(J) - h*d*J\n", - " else\n", - " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", - " # M * dy/dt = F(t,y)\n", - " # we can simply replace eye(J) by M in the following expression\n", - " # (see Sec. 5 in [SR97])\n", - "\n", - " W = lufact( eye(J) - h*d*J )\n", - " end\n", - "\n", - " # approximate time-derivative of F\n", - " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", - "\n", - " # modified Rosenbrock formula\n", - " k1 = W\\(F0 + T)\n", - " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", - " k2 = W\\(F1 - k1) + k1\n", - " ynew = y + h*k2\n", - " F2 = F(t + h, ynew)\n", - " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", - "\n", - " \n", - " # Component-wise error check similar to MATLAB ode23s\n", - " \n", - " threshold = abstol/reltol\n", - " \n", - " if normControl\n", - " err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate\n", - " else\n", - " err = (abs(h)/6)*(norm(k1 - 2*k2 + k3)/max(max(norm(y),norm(ynew)),threshold)) # error estimate\n", - " end\n", - " \n", - " # check if new solution is acceptable\n", - " if err <= reltol\n", - "\n", - " if points==:specified || points==:all\n", - " # only points in tspan are requested\n", - " # -> find relevant points in (t,t+h]\n", - " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", - " # rescale to (0,1]\n", - " s = (toi-t)/h\n", - "\n", - " # use interpolation formula to get solutions at t=toi\n", - " push!(tout, toi)\n", - " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", - " end\n", - " end\n", - " if (points==:all) && (tout[end]!=t+h)\n", - " # add the intermediate points\n", - " push!(tout, t + h)\n", - " push!(yout, ynew)\n", - " end\n", - "\n", - " # update solution\n", - " t = t + h\n", - " y = ynew\n", - "\n", - " F0 = F2 # use FSAL property\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - " # for new solution\n", - " end\n", - "\n", - " # update of the step size\n", - " h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) )\n", - " end\n", - "\n", - " return tout, yout\n", - "end\n", - "\n", - "\n", - "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", - "# with provided coefficients.\n", - "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", - "\n", - " if typeof(jacobian) == Function\n", - " G = jacobian\n", - " else\n", - " G = (t, x)->fdjacobian(F, x, t)\n", - " end\n", - "\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " solstep = 1\n", - " while solstep < length(tspan)\n", - " ts = tspan[solstep]\n", - " hs = h[solstep]\n", - " xs = x[solstep]\n", - " dFdx = G(ts, xs)\n", - " # FIXME\n", - " if size(dFdx,1) == 1\n", - " jac = 1/gamma/hs - dFdx[1]\n", - " else\n", - " jac = eye(dFdx)/gamma/hs - dFdx\n", - " end\n", - "\n", - " g = Array(typeof(x0), size(a,1))\n", - " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", - " x[solstep+1] = x[solstep] + b[1]*g[1]\n", - "\n", - " for i = 2:size(a,1)\n", - " dx = zero(x0)\n", - " dF = zero(x0/hs)\n", - " for j = 1:i-1\n", - " dx += a[i,j]*g[j]\n", - " dF += c[i,j]*g[j]\n", - " end\n", - " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", - " x[solstep+1] += b[i]*g[i]\n", - " end\n", - " solstep += 1\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "\n", - "# Kaps-Rentrop coefficients\n", - "const kr4_coefficients = (0.231,\n", - " [0 0 0 0\n", - " 2 0 0 0\n", - " 4.452470820736 4.16352878860 0 0\n", - " 4.452470820736 4.16352878860 0 0],\n", - " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", - " [ 0 0 0 0\n", - " -5.07167533877 0 0 0\n", - " 6.02015272865 0.1597500684673 0 0\n", - " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", - "\n", - "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Shampine coefficients\n", - "const s4_coefficients = (0.5,\n", - " [ 0 0 0 0\n", - " 2 0 0 0\n", - " 48/25 6/25 0 0\n", - " 48/25 6/25 0 0],\n", - " [19/9 1/2 25/108 125/108],\n", - " [ 0 0 0 0\n", - " -8 0 0 0\n", - " 372/25 12/5 0 0\n", - " -112/125 -54/125 -2/5 0],)\n", - "\n", - "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Use Shampine coefficients by default (matching Numerical Recipes)\n", - "const ode4s = ode4s_s\n", - "\n", - "const ms_coefficients4 = [ 1 0 0 0\n", - " -1/2 3/2 0 0\n", - " 5/12 -4/3 23/12 0\n", - " -9/24 37/24 -59/24 55/24]\n", - "\n", - "\n", - "#end # module ODE\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "tol = 1e-2\n", - "function norm1(x,p=1)\n", - " return norm(x,p)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;])\n", - "ys = hcat(y...).'\n", - "maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->y, 1., [1:-.001:0;])\n", - "maximum(abs(y-e.^(t-1))) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->y, 1., [0:.001:1;])\n", - "maximum(abs(y-e.^t)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->2t, 0., [0:.001:1;])\n", - "maximum(abs(y-t.^2)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->6.0, 0., [0:.1:1;];norm=norm1)\n", - "maximum(abs(y-6t)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "methods(norm)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 0.4.2", - "language": "julia", - "name": "julia-0.4" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "0.4.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} From 444468e2070b79fd2b3d540c4bcc61a892ddbf62 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 17:24:43 +0200 Subject: [PATCH 04/14] Updated the files with comments --- src/ODE.jl | 19 +- src/Untitled.ipynb | 582 --------------------------------------------- src/runge_kutta.jl | 10 +- 3 files changed, 18 insertions(+), 593 deletions(-) delete mode 100644 src/Untitled.ipynb diff --git a/src/ODE.jl b/src/ODE.jl index 52374c8e9..caf4e851f 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -67,7 +67,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) tdir = sign(tend-t0) tdir==0 && error("Zero time span") tau = max(reltol.*abs(x0), abstol) - d0 = norm(x0./tau, Inf) + d0 = norm(x0,./tau Inf) f0 = F(t0, x0) d1 = norm(f0./tau, Inf) if d0 < 1e-5 || d1 < 1e-5 @@ -247,7 +247,14 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # constants const d = 1/(2 + sqrt(2)) const e32 = 6 + sqrt(2) + + + # initialization + t = tspan[1] + # Component-wise Tolerance check similar to MATLAB ode23s + # Support for component-wise absolute tolerance only. + # Relative tolerance should be a scalar. @assert length(reltol) == 1 "Relative tolerance must be a scalar" @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" @@ -255,9 +262,6 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, if length(abstol) == 1 && length(y0) != 1 abstol = abstol*ones(y0); end - - # initialization - t = tspan[1] tfinal = tspan[end] @@ -307,12 +311,9 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, F2 = F(t + h, ynew) k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) - - # Component-wise error check similar to MATLAB ode23s - - threshold = abstol/reltol + threshold = abstol/reltol # error threshold err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate - + # check if new solution is acceptable if err <= reltol diff --git a/src/Untitled.ipynb b/src/Untitled.ipynb deleted file mode 100644 index 8da022fb6..000000000 --- a/src/Untitled.ipynb +++ /dev/null @@ -1,582 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "#isdefined(Base, :__precompile__) && __precompile__()\n", - "# Ordinary Differential Equation Solvers\n", - "\n", - "#module ODE\n", - "\n", - "using Polynomials\n", - "using Compat\n", - "\n", - "## minimal function export list\n", - "# adaptive non-stiff:\n", - "export ode23, ode45, ode78\n", - "# non-adaptive non-stiff:\n", - "export ode4, ode4ms\n", - "# adaptive stiff:\n", - "export ode23s\n", - "# non-adaptive stiff:\n", - "export ode4s\n", - "\n", - "## complete function export list: see runtests.jl\n", - "\n", - "###############################################################################\n", - "## Coefficient Tableaus\n", - "###############################################################################\n", - "\n", - "# Butcher Tableaus, or more generally coefficient tables\n", - "# see Hairer & Wanner 1992, p. 134, 166\n", - "\n", - "abstract Tableau{Name, S, T<:Real}\n", - "# Name is the name of the tableau/method (a symbol)\n", - "# S is the number of stages (an int)\n", - "# T is the type of the coefficients\n", - "#\n", - "# TODO: have a type parameter which specifies adaptive vs non-adaptive\n", - "#\n", - "# For all types of tableaus it assumes fields:\n", - "# order::(Int...) # order of the method(s)\n", - "#\n", - "# For Runge-Kutta methods it assumes fields:\n", - "# a::Matrix{T} # SxS matrix\n", - "# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive)\n", - "# c::Vector{T} # S\n", - "#\n", - "# For a tableau:\n", - "# c1 | a_11 .... a_1s\n", - "# . | a_21 . .\n", - "# . | a_31 . .\n", - "# . | .... . .\n", - "# c_s | a_s1 ....... a_ss\n", - "# -----+--------------------\n", - "# | b_1 ... b_s this is the one used for stepping\n", - "# | b'_1 ... b'_s this is the one used for error-checking\n", - "\n", - "Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T\n", - "order(b::Tableau) = b.order\n", - "# Subtypes need to define a convert method to convert to a different\n", - "# eltype with signature:\n", - "Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error(\"Define convert method for concrete Tableau types\")\n", - "\n", - "###############################################################################\n", - "## HELPER FUNCTIONS\n", - "###############################################################################\n", - "\n", - "# estimator for initial step based on book\n", - "# \"Solving Ordinary Differential Equations I\" by Hairer et al., p.169\n", - "function hinit(F, x0, t0, tend, p, reltol, abstol)\n", - " # Returns first step, direction of integration and F evaluated at t0\n", - " tdir = sign(tend-t0)\n", - " tdir==0 && error(\"Zero time span\")\n", - " tau = max(reltol.*abs(x0), abstol)\n", - " d0 = norm(x0./tau, Inf)\n", - " f0 = F(t0, x0)\n", - " d1 = norm(f0./tau, Inf)\n", - " if d0 < 1e-5 || d1 < 1e-5\n", - " h0 = 1e-6\n", - " else\n", - " h0 = 0.01*(d0/d1)\n", - " end\n", - " # perform Euler step\n", - " x1 = x0 + tdir*h0*f0\n", - " f1 = F(t0 + tdir*h0, x1)\n", - " # estimate second derivative\n", - " d2 = norm((f1 - f0)./tau, Inf)/h0\n", - " if max(d1, d2) <= 1e-15\n", - " h1 = max(1e-6, 1e-3*h0)\n", - " else\n", - " pow = -(2. + log10(max(d1, d2)))/(p + 1.)\n", - " h1 = 10.^pow\n", - " end\n", - " return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0\n", - "end\n", - "\n", - "# isoutofdomain takes the state and returns true if state is outside\n", - "# of the allowed domain. Used in adaptive step-control.\n", - "isoutofdomain(x) = isnan(x)\n", - "\n", - "function make_consistent_types(fn, y0, tspan, btab::Tableau)\n", - " # There are a few types involved in a call to a ODE solver which\n", - " # somehow need to be consistent:\n", - " #\n", - " # Et = eltype(tspan)\n", - " # Ey = eltype(y0)\n", - " # Ef = eltype(Tf)\n", - " #\n", - " # There are also the types of the containers, but they are not\n", - " # needed as `similar` is used to make containers.\n", - " # Tt = typeof(tspan)\n", - " # Ty = typeof(y0) # note, this can be a scalar\n", - " # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar\n", - " #\n", - " # Returns\n", - " # - Et: eltype of time, needs to be a real \"continuous\" type, at\n", - " # the moment a AbstractFloat\n", - " # - Eyf: suitable eltype of y and f(t,y)\n", - " # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))\n", - " # - Ty: container type of y0\n", - " # - btab: tableau with entries converted to Et\n", - "\n", - " # Needed interface:\n", - " # On components: /, -\n", - " # On container: eltype, promote_type\n", - " # On time container: eltype\n", - "\n", - " Ty = typeof(y0)\n", - " Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))\n", - "\n", - " Et = eltype(tspan)\n", - " @assert Et<:Real\n", - " if !(Et<:AbstractFloat)\n", - " Et = promote_type(Et, Float64)\n", - " end\n", - "\n", - " # if all are Floats, make them the same\n", - " if Et<:AbstractFloat && Eyf<:AbstractFloat\n", - " Et = promote_type(Et, Eyf)\n", - " Eyf = Et\n", - " end\n", - "\n", - " !isleaftype(Et) && warn(\"The eltype(tspan) is not a concrete type! Change type of tspan for better performance.\")\n", - " !isleaftype(Eyf) && warn(\"The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.\")\n", - "\n", - " btab_ = convert(Et, btab)\n", - " return Et, Eyf, Ty, btab_\n", - "end\n", - "\n", - "###############################################################################\n", - "## NON-STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "include(\"runge_kutta.jl\")\n", - "\n", - "# ODE_MS Fixed-step, fixed-order multi-step numerical method\n", - "# with Adams-Bashforth-Moulton coefficients\n", - "function ode_ms(F, x0, tspan, order::Integer)\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " if 1 <= order <= 4\n", - " b = ms_coefficients4\n", - " else\n", - " b = zeros(order, order)\n", - " b[1:4, 1:4] = ms_coefficients4\n", - " for s = 5:order\n", - " for j = 0:(s - 1)\n", - " # Assign in correct order for multiplication below\n", - " # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)\n", - " p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))\n", - " b[s, s - j] = ((-1)^j / factorial(j)\n", - " / factorial(s - 1 - j) * polyval(p_int, 1))\n", - " end\n", - " end\n", - " end\n", - "\n", - " # TODO: use a better data structure here (should be an order-element circ buffer)\n", - " xdot = similar(x)\n", - " for i = 1:length(tspan)-1\n", - " # Need to run the first several steps at reduced order\n", - " steporder = min(i, order)\n", - " xdot[i] = F(tspan[i], x[i])\n", - "\n", - " x[i+1] = x[i]\n", - " for j = 1:steporder\n", - " x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]\n", - " end\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "# Use order 4 by default\n", - "ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)\n", - "ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)\n", - "\n", - "###############################################################################\n", - "## STIFF SOLVERS\n", - "###############################################################################\n", - "\n", - "# Crude forward finite differences estimator of Jacobian as fallback\n", - "\n", - "# FIXME: This doesn't really work if x is anything but a Vector or a scalar\n", - "function fdjacobian(F, x::Number, t)\n", - " ftx = F(t, x)\n", - "\n", - " # The 100 below is heuristic\n", - " dx = (x .+ (x==0))./100\n", - " dFdx = (F(t,x+dx)-ftx)./dx\n", - "\n", - " return dFdx\n", - "end\n", - "\n", - "function fdjacobian(F, x, t)\n", - " ftx = F(t, x)\n", - " lx = max(length(x),1)\n", - " dFdx = zeros(eltype(x), lx, lx)\n", - " for j = 1:lx\n", - " # The 100 below is heuristic\n", - " dx = zeros(eltype(x), lx)\n", - " dx[j] = (x[j] .+ (x[j]==0))./100\n", - " dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]\n", - " end\n", - " return dFdx\n", - "end\n", - "\n", - "# ODE23S Solve stiff systems based on a modified Rosenbrock triple\n", - "# (also used by MATLAB's ODE23s); see Sec. 4.1 in\n", - "#\n", - "# [SR97] L.F. Shampine and M.W. Reichelt: \"The MATLAB ODE Suite,\" SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22\n", - "#\n", - "# supports keywords: points = :all | :specified (using dense output)\n", - "# jacobian = G(t,y)::Function | nothing (FD)\n", - "function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,\n", - " jacobian=nothing,\n", - " points=:all,\n", - " norm=Base.norm,\n", - " minstep=abs(tspan[end] - tspan[1])/1e18,\n", - " maxstep=abs(tspan[end] - tspan[1])/2.5,\n", - " initstep=0.)\n", - "\n", - "\n", - " # select method for computing the Jacobian\n", - " if typeof(jacobian) == Function\n", - " jac = jacobian\n", - " else\n", - " # fallback finite-difference\n", - " jac = (t, y)->fdjacobian(F, y, t)\n", - " end\n", - "\n", - " # constants\n", - " const d = 1/(2 + sqrt(2))\n", - " const e32 = 6 + sqrt(2)\n", - "\n", - " # Check if norm provided by the user\n", - " normControl = (norm == Base.norm)\n", - " \n", - " @assert length(reltol) == 1 \"Relative tolerance must be a scalar\"\n", - " @assert length(abstol) == 1 || length(abstol) == length(y0) \"Dimension of Absolute tolerance does not match the dimension of the problem\"\n", - " if !normControl\n", - " @assert length(abstol) == 1 \"Norm control supports only scalar Absolute tolerance\"\n", - " end\n", - " \n", - " # Broadcast the abstol to a vector\n", - " if length(abstol) == 1 && length(y0) != 1 && normControl\n", - " abstol = abstol*ones(y0);\n", - " end\n", - " \n", - " # initialization\n", - " t = tspan[1]\n", - "\n", - " tfinal = tspan[end]\n", - "\n", - " h = initstep\n", - " if h == 0.\n", - " # initial guess at a step size\n", - " h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol)\n", - " else\n", - " tdir = sign(tfinal - t)\n", - " F0 = F(t,y0)\n", - " end\n", - " h = tdir * min(abs(h), maxstep)\n", - "\n", - " y = y0\n", - " tout = Array(typeof(t), 1)\n", - " tout[1] = t # first output time\n", - " yout = Array(typeof(y0), 1)\n", - " yout[1] = deepcopy(y) # first output solution\n", - "\n", - "\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - "\n", - " while abs(t - tfinal) > 0 && minstep < abs(h)\n", - " if abs(t-tfinal) < abs(h)\n", - " h = tfinal - t\n", - " end\n", - "\n", - " if size(J,1) == 1\n", - " W = one(J) - h*d*J\n", - " else\n", - " # note: if there is a mass matrix M on the lhs of the ODE, i.e.,\n", - " # M * dy/dt = F(t,y)\n", - " # we can simply replace eye(J) by M in the following expression\n", - " # (see Sec. 5 in [SR97])\n", - "\n", - " W = lufact( eye(J) - h*d*J )\n", - " end\n", - "\n", - " # approximate time-derivative of F\n", - " T = h*d*(F(t + h/100, y) - F0)/(h/100)\n", - "\n", - " # modified Rosenbrock formula\n", - " k1 = W\\(F0 + T)\n", - " F1 = F(t + 0.5*h, y + 0.5*h*k1)\n", - " k2 = W\\(F1 - k1) + k1\n", - " ynew = y + h*k2\n", - " F2 = F(t + h, ynew)\n", - " k3 = W\\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )\n", - "\n", - " \n", - " # Component-wise error check similar to MATLAB ode23s\n", - " \n", - " threshold = abstol/reltol\n", - " \n", - " if normControl\n", - " err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate\n", - " else\n", - " err = (abs(h)/6)*(norm(k1 - 2*k2 + k3)/max(max(norm(y),norm(ynew)),threshold)) # error estimate\n", - " end\n", - " \n", - " # check if new solution is acceptable\n", - " if err <= reltol\n", - "\n", - " if points==:specified || points==:all\n", - " # only points in tspan are requested\n", - " # -> find relevant points in (t,t+h]\n", - " for toi in tspan[(tspan.>t) & (tspan.<=t+h)]\n", - " # rescale to (0,1]\n", - " s = (toi-t)/h\n", - "\n", - " # use interpolation formula to get solutions at t=toi\n", - " push!(tout, toi)\n", - " push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))\n", - " end\n", - " end\n", - " if (points==:all) && (tout[end]!=t+h)\n", - " # add the intermediate points\n", - " push!(tout, t + h)\n", - " push!(yout, ynew)\n", - " end\n", - "\n", - " # update solution\n", - " t = t + h\n", - " y = ynew\n", - "\n", - " F0 = F2 # use FSAL property\n", - " J = jac(t,y) # get Jacobian of F wrt y\n", - " # for new solution\n", - " end\n", - "\n", - " # update of the step size\n", - " h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) )\n", - " end\n", - "\n", - " return tout, yout\n", - "end\n", - "\n", - "\n", - "#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method\n", - "# with provided coefficients.\n", - "function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)\n", - "\n", - " if typeof(jacobian) == Function\n", - " G = jacobian\n", - " else\n", - " G = (t, x)->fdjacobian(F, x, t)\n", - " end\n", - "\n", - " h = diff(tspan)\n", - " x = Array(typeof(x0), length(tspan))\n", - " x[1] = x0\n", - "\n", - " solstep = 1\n", - " while solstep < length(tspan)\n", - " ts = tspan[solstep]\n", - " hs = h[solstep]\n", - " xs = x[solstep]\n", - " dFdx = G(ts, xs)\n", - " # FIXME\n", - " if size(dFdx,1) == 1\n", - " jac = 1/gamma/hs - dFdx[1]\n", - " else\n", - " jac = eye(dFdx)/gamma/hs - dFdx\n", - " end\n", - "\n", - " g = Array(typeof(x0), size(a,1))\n", - " g[1] = (jac \\ F(ts + b[1]*hs, xs))\n", - " x[solstep+1] = x[solstep] + b[1]*g[1]\n", - "\n", - " for i = 2:size(a,1)\n", - " dx = zero(x0)\n", - " dF = zero(x0/hs)\n", - " for j = 1:i-1\n", - " dx += a[i,j]*g[j]\n", - " dF += c[i,j]*g[j]\n", - " end\n", - " g[i] = (jac \\ (F(ts + b[i]*hs, xs + dx) + dF/hs))\n", - " x[solstep+1] += b[i]*g[i]\n", - " end\n", - " solstep += 1\n", - " end\n", - " return vcat(tspan), x\n", - "end\n", - "\n", - "\n", - "# Kaps-Rentrop coefficients\n", - "const kr4_coefficients = (0.231,\n", - " [0 0 0 0\n", - " 2 0 0 0\n", - " 4.452470820736 4.16352878860 0 0\n", - " 4.452470820736 4.16352878860 0 0],\n", - " [3.95750374663 4.62489238836 0.617477263873 1.28261294568],\n", - " [ 0 0 0 0\n", - " -5.07167533877 0 0 0\n", - " 6.02015272865 0.1597500684673 0 0\n", - " -1.856343618677 -8.50538085819 -2.08407513602 0],)\n", - "\n", - "ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Shampine coefficients\n", - "const s4_coefficients = (0.5,\n", - " [ 0 0 0 0\n", - " 2 0 0 0\n", - " 48/25 6/25 0 0\n", - " 48/25 6/25 0 0],\n", - " [19/9 1/2 25/108 125/108],\n", - " [ 0 0 0 0\n", - " -8 0 0 0\n", - " 372/25 12/5 0 0\n", - " -112/125 -54/125 -2/5 0],)\n", - "\n", - "ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian)\n", - "\n", - "# Use Shampine coefficients by default (matching Numerical Recipes)\n", - "const ode4s = ode4s_s\n", - "\n", - "const ms_coefficients4 = [ 1 0 0 0\n", - " -1/2 3/2 0 0\n", - " 5/12 -4/3 23/12 0\n", - " -9/24 37/24 -59/24 55/24]\n", - "\n", - "\n", - "#end # module ODE\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "tol = 1e-2\n", - "function norm1(x,p=1)\n", - " return norm(x,p)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;])\n", - "ys = hcat(y...).'\n", - "maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->y, 1., [1:-.001:0;])\n", - "maximum(abs(y-e.^(t-1))) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->y, 1., [0:.001:1;])\n", - "maximum(abs(y-e.^t)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->2t, 0., [0:.001:1;])\n", - "maximum(abs(y-t.^2)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "t,y=ode23s((t,y)->6.0, 0., [0:.1:1;];norm=norm1)\n", - "maximum(abs(y-6t)) < tol" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "methods(norm)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 0.4.2", - "language": "julia", - "name": "julia-0.4" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "0.4.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 5f1a994a6..649a6e137 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -251,12 +251,18 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici tstart = tspan[1] tend = tspan[end] - if length(reltol) == 1 + # Component-wise reltol and abstol as per Solving Ordinary Differential Equations I by Hairer and Wanner + @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" + + # Broadcast the tolerances if scalars are provided + if length(reltol) == 1 && length(y0) != 1 reltol = reltol*ones(y0); end - if length(abstol) == 1 + if length(abstol) == 1 && length(y0) != 1 abstol = abstol*ones(y0); end + # work arrays: y = similar(y0, Eyf, dof) # y at time t From 23349f38637a29fc7ebbeec52e0f1f7980784529 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 20:38:07 +0200 Subject: [PATCH 05/14] Fixed the error in line 70 in ODE.jl --- src/ODE.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ODE.jl b/src/ODE.jl index caf4e851f..64d1db18d 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -67,7 +67,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) tdir = sign(tend-t0) tdir==0 && error("Zero time span") tau = max(reltol.*abs(x0), abstol) - d0 = norm(x0,./tau Inf) + d0 = norm(x0./tau, Inf) f0 = F(t0, x0) d1 = norm(f0./tau, Inf) if d0 < 1e-5 || d1 < 1e-5 From 740b05eea531f053ab4781ac1dfb9b2f5b5aa40a Mon Sep 17 00:00:00 2001 From: sonVishal Date: Thu, 31 Mar 2016 20:41:32 +0200 Subject: [PATCH 06/14] Added Travis file as required --- .travis.yml | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 .travis.yml diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 000000000..1fc28491b --- /dev/null +++ b/.travis.yml @@ -0,0 +1,17 @@ +language: julia +os: + - osx + - linux +julia: + - 0.3 + - release + - nightly +git: + depth: 999999 +notifications: + email: false +script: + - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi + - julia -e 'Pkg.clone(pwd()); Pkg.build("ODE"); Pkg.test("ODE"; coverage=true)'; +after_success: + - julia -e 'cd(Pkg.dir("ODE")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'; From 58bd9cacf409d9375046f50f86a66b412f8b09d4 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 08:11:51 +0200 Subject: [PATCH 07/14] Corrected the runtest.jl Robertson test case and removed the comments. Corrected reltol::Number for ode23s --- src/ODE.jl | 5 ++--- test/runtests.jl | 8 ++------ 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 64d1db18d..7c3fab792 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -227,7 +227,7 @@ end # # supports keywords: points = :all | :specified (using dense output) # jacobian = G(t,y)::Function | nothing (FD) -function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, +function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, jacobian=nothing, points=:all, norm=Base.norm, @@ -255,7 +255,6 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # Component-wise Tolerance check similar to MATLAB ode23s # Support for component-wise absolute tolerance only. # Relative tolerance should be a scalar. - @assert length(reltol) == 1 "Relative tolerance must be a scalar" @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" # Broadcast the abstol to a vector @@ -312,7 +311,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) threshold = abstol/reltol # error threshold - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # scaled error estimate # check if new solution is acceptable if err <= reltol diff --git a/test/runtests.jl b/test/runtests.jl index c2c211178..b66073abc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,10 +65,6 @@ end # rober testcase from http://www.unige.ch/~hairer/testset/testset.html -# Changing the test value. -# Even the Matlab ode23s does not give the error to be less than 2e-10. -# The value comes out to be approximately 2.042354e-10. -# Hence changing the test value to 2.1e-10 let println("ROBER test case") function f(t, y) @@ -79,13 +75,13 @@ let ydot end t = [0., 1e11] - t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, + t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-9, reltol=1e-9, maxstep=1e11/10, minstep=1e11/1e18) refsol = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050] # reference solution at tspan[2] - @test norm(refsol-y[end], Inf) < 2.1e-10 + @test norm(refsol-y[end], Inf) < 2e-10 end include("interface-tests.jl") From 936931a3b287f9464357eda9e3d86d64535efe3a Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 08:25:56 +0200 Subject: [PATCH 08/14] Was using Inf norm for err. Now using only norm() since it is provided by the user and default value is provided. --- src/ODE.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ODE.jl b/src/ODE.jl index 7c3fab792..6087f48c4 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -311,7 +311,7 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) threshold = abstol/reltol # error threshold - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # scaled error estimate + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold)) # scaled error estimate # check if new solution is acceptable if err <= reltol From 9c89efb8193f5b637a2e436134c99c39840ddd41 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 09:10:24 +0200 Subject: [PATCH 09/14] Component-wise reltol for ode23s. In case of scalar reltol "and" abstol, restore old behaviour. --- src/ODE.jl | 36 ++++++++++++++++++------------------ src/runge_kutta.jl | 33 ++++++++++++++++++--------------- 2 files changed, 36 insertions(+), 33 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 6087f48c4..fd7002146 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,14 +62,14 @@ Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert met # estimator for initial step based on book # "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0, tend, p, reltol, abstol) +function hinit(F, x0, t0, tend, p, reltol, abstol, norm) # Returns first step, direction of integration and F evaluated at t0 tdir = sign(tend-t0) tdir==0 && error("Zero time span") tau = max(reltol.*abs(x0), abstol) - d0 = norm(x0./tau, Inf) + d0 = norm(x0./tau) f0 = F(t0, x0) - d1 = norm(f0./tau, Inf) + d1 = norm(f0./tau) if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 else @@ -79,7 +79,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) x1 = x0 + tdir*h0*f0 f1 = F(t0 + tdir*h0, x1) # estimate second derivative - d2 = norm((f1 - f0)./tau, Inf)/h0 + d2 = norm((f1 - f0)./tau)/h0 if max(d1, d2) <= 1e-15 h1 = max(1e-6, 1e-3*h0) else @@ -227,7 +227,7 @@ end # # supports keywords: points = :all | :specified (using dense output) # jacobian = G(t,y)::Function | nothing (FD) -function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, +function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, jacobian=nothing, points=:all, norm=Base.norm, @@ -252,22 +252,16 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, # initialization t = tspan[1] - # Component-wise Tolerance check similar to MATLAB ode23s - # Support for component-wise absolute tolerance only. - # Relative tolerance should be a scalar. + # Component-wise Tolerance check similar to the non-stiff solvers @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" - - # Broadcast the abstol to a vector - if length(abstol) == 1 && length(y0) != 1 - abstol = abstol*ones(y0); - end + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" tfinal = tspan[end] h = initstep if h == 0. # initial guess at a step size - h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) + h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol, norm) else tdir = sign(tfinal - t) F0 = F(t,y0) @@ -310,11 +304,17 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, F2 = F(t + h, ynew) k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) - threshold = abstol/reltol # error threshold - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold)) # scaled error estimate + # If reltol and abstol are vectors + # perform component-wise computations for error + if length(abstol) != 1 && length(reltol) != 1 + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate + else + # else restore old behvaiour + err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate + end # check if new solution is acceptable - if err <= reltol + if err <= 1 if points==:specified || points==:all # only points in tspan are requested @@ -344,7 +344,7 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, end # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) ) + h = tdir*min( maxstep, abs(h)*0.8*(1/err)^(1/3) ) end return tout, yout diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 649a6e137..651880c08 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -253,16 +253,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici # Component-wise reltol and abstol as per Solving Ordinary Differential Equations I by Hairer and Wanner @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" - @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" - - # Broadcast the tolerances if scalars are provided - if length(reltol) == 1 && length(y0) != 1 - reltol = reltol*ones(y0); - end - if length(abstol) == 1 && length(y0) != 1 - abstol = abstol*ones(y0); - end - + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" # work arrays: y = similar(y0, Eyf, dof) # y at time t @@ -289,7 +280,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici error("Unrecognized option points==$points") end # Time - dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 + dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol, norm) # sets ks[1]=f0 if initstep!=0 dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") end @@ -422,10 +413,22 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, facmin = 1./facmax # maximal step size decrease. ? # in-place calculate xerr./tol - for d=1:dof - # if outside of domain (usually NaN) then make step size smaller by maximum - isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 + + # If reltol and abstol are vectors + # perform component-wise computations + if length(abstol) != 1 && length(reltol) != 1 + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan + xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 + end + else + # else restore old behaviour + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan + xerr[d] = xerr[d]/(abstol + max(abs(x0[d]), abs(xtrial[d]))*reltol) # Eq 4.10 + end end err = norm(xerr, 2) # Eq. 4.11 newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified From 84d01ad6482cd84c334e2e842ee9f830d3c9a530 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 09:10:24 +0200 Subject: [PATCH 10/14] Component-wise reltol for ode23s. In case of scalar reltol "and" abstol, restore old behaviour. Also, passing "norm" to hinit so that the user defined norm can be used for the initial step size computation. --- src/ODE.jl | 36 ++++++++++++++++++------------------ src/runge_kutta.jl | 33 ++++++++++++++++++--------------- 2 files changed, 36 insertions(+), 33 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 6087f48c4..fd7002146 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,14 +62,14 @@ Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert met # estimator for initial step based on book # "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0, tend, p, reltol, abstol) +function hinit(F, x0, t0, tend, p, reltol, abstol, norm) # Returns first step, direction of integration and F evaluated at t0 tdir = sign(tend-t0) tdir==0 && error("Zero time span") tau = max(reltol.*abs(x0), abstol) - d0 = norm(x0./tau, Inf) + d0 = norm(x0./tau) f0 = F(t0, x0) - d1 = norm(f0./tau, Inf) + d1 = norm(f0./tau) if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 else @@ -79,7 +79,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) x1 = x0 + tdir*h0*f0 f1 = F(t0 + tdir*h0, x1) # estimate second derivative - d2 = norm((f1 - f0)./tau, Inf)/h0 + d2 = norm((f1 - f0)./tau)/h0 if max(d1, d2) <= 1e-15 h1 = max(1e-6, 1e-3*h0) else @@ -227,7 +227,7 @@ end # # supports keywords: points = :all | :specified (using dense output) # jacobian = G(t,y)::Function | nothing (FD) -function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, +function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, jacobian=nothing, points=:all, norm=Base.norm, @@ -252,22 +252,16 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, # initialization t = tspan[1] - # Component-wise Tolerance check similar to MATLAB ode23s - # Support for component-wise absolute tolerance only. - # Relative tolerance should be a scalar. + # Component-wise Tolerance check similar to the non-stiff solvers @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" - - # Broadcast the abstol to a vector - if length(abstol) == 1 && length(y0) != 1 - abstol = abstol*ones(y0); - end + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" tfinal = tspan[end] h = initstep if h == 0. # initial guess at a step size - h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) + h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol, norm) else tdir = sign(tfinal - t) F0 = F(t,y0) @@ -310,11 +304,17 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, F2 = F(t + h, ynew) k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) - threshold = abstol/reltol # error threshold - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold)) # scaled error estimate + # If reltol and abstol are vectors + # perform component-wise computations for error + if length(abstol) != 1 && length(reltol) != 1 + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate + else + # else restore old behvaiour + err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate + end # check if new solution is acceptable - if err <= reltol + if err <= 1 if points==:specified || points==:all # only points in tspan are requested @@ -344,7 +344,7 @@ function ode23s(F, y0, tspan; reltol::Number = 1.0e-5, abstol = 1.0e-8, end # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) ) + h = tdir*min( maxstep, abs(h)*0.8*(1/err)^(1/3) ) end return tout, yout diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 649a6e137..651880c08 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -253,16 +253,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici # Component-wise reltol and abstol as per Solving Ordinary Differential Equations I by Hairer and Wanner @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" - @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" - - # Broadcast the tolerances if scalars are provided - if length(reltol) == 1 && length(y0) != 1 - reltol = reltol*ones(y0); - end - if length(abstol) == 1 && length(y0) != 1 - abstol = abstol*ones(y0); - end - + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" # work arrays: y = similar(y0, Eyf, dof) # y at time t @@ -289,7 +280,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici error("Unrecognized option points==$points") end # Time - dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 + dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol, norm) # sets ks[1]=f0 if initstep!=0 dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") end @@ -422,10 +413,22 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, facmin = 1./facmax # maximal step size decrease. ? # in-place calculate xerr./tol - for d=1:dof - # if outside of domain (usually NaN) then make step size smaller by maximum - isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 + + # If reltol and abstol are vectors + # perform component-wise computations + if length(abstol) != 1 && length(reltol) != 1 + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan + xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 + end + else + # else restore old behaviour + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan + xerr[d] = xerr[d]/(abstol + max(abs(x0[d]), abs(xtrial[d]))*reltol) # Eq 4.10 + end end err = norm(xerr, 2) # Eq. 4.11 newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified From 76a8e05dc0c0a9981d5a2a3f190f8a7babf95c30 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 10:01:19 +0200 Subject: [PATCH 11/14] Had written the if statements the other way around. Should execute the old behaviour only if both abstol and reltol are scalars. If either of them are vectors, perform component-wise computations. --- src/ODE.jl | 10 +++++----- src/runge_kutta.jl | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index fd7002146..d66b734a9 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -305,12 +305,12 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) # If reltol and abstol are vectors - # perform component-wise computations for error - if length(abstol) != 1 && length(reltol) != 1 - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate - else - # else restore old behvaiour + # restore old behvaiour + if length(abstol) == 1 && length(reltol) == 1 err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate + else + # else perform component-wise computations for error + err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate end # check if new solution is acceptable diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index 651880c08..a43933404 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -415,19 +415,19 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, # in-place calculate xerr./tol # If reltol and abstol are vectors - # perform component-wise computations - if length(abstol) != 1 && length(reltol) != 1 + # restore old behaviour + if length(abstol) == 1 && length(reltol) == 1 for d=1:dof # if outside of domain (usually NaN) then make step size smaller by maximum isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 + xerr[d] = xerr[d]/(abstol + max(abs(x0[d]), abs(xtrial[d]))*reltol) # Eq 4.10 end else - # else restore old behaviour + # else perform component-wise computations for d=1:dof # if outside of domain (usually NaN) then make step size smaller by maximum isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol + max(abs(x0[d]), abs(xtrial[d]))*reltol) # Eq 4.10 + xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 end end err = norm(xerr, 2) # Eq. 4.11 From ea10a9d51c6f04e0f62eee6c9ad36d12ca09f111 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Fri, 1 Apr 2016 17:58:18 +0200 Subject: [PATCH 12/14] Implemented getScaledError as a multiple dispatch function for both stiff and non-stiff solver This avoids the 'if' construct --- src/ODE.jl | 46 +++++++++++++++++++++++++++++++++++++++++----- src/runge_kutta.jl | 18 ++---------------- 2 files changed, 43 insertions(+), 21 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index d66b734a9..93f5745fb 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -220,6 +220,40 @@ function fdjacobian(F, x, t) return dFdx end +## Function to get the Scaled error given the error estimate +# The following 3 functions are for the stiff solver +# If both reltol and abstol are scalars then restore previous behaviour +function getScaledError(y,ynew,errEstimate,h,reltol::Number,abstol::Number,norm) + return (abs(h)/6)*norm(errEstimate)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate +end +# If atleast one of them is a vector then perform component-wise computations +function getScaledError(y,ynew,errEstimate,h,reltol::Number,abstol::Vector,norm) + return (abs(h)/6)*norm((errEstimate)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate +end +function getScaledError(y,ynew,errEstimate,h,reltol::Vector,abstol,norm) + return (abs(h)/6)*norm((errEstimate)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate +end +# The following 2 functions are for the non-stiff solvers since it requires dof check +# As per the old code, the error estimate is rewritten with the scaled error +# Estimates the error and a new step size following Hairer & +# Wanner 1992, p167 (with some modifications) +function getScaledError!(y,ynew,errEstimate,reltol::Number,abstol::Number,norm,dof) + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(ynew[d]) && return true + errEstimate[d] = errEstimate[d]/(abstol + max(abs(y[d]), abs(ynew[d]))*reltol) # Eq 4.10 + end + return false +end +function getScaledError!(y,ynew,errEstimate,reltol::Vector,abstol::Vector,norm,dof) + for d=1:dof + # if outside of domain (usually NaN) then break and return + isoutofdomain(ynew[d]) && return true + errEstimate[d] = errEstimate[d]/(abstol[d] + max(abs(y[d]), abs(ynew[d]))*reltol[d]) # Eq 4.10 + end + return false +end + # ODE23S Solve stiff systems based on a modified Rosenbrock triple # (also used by MATLAB's ODE23s); see Sec. 4.1 in # @@ -306,12 +340,14 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # If reltol and abstol are vectors # restore old behvaiour - if length(abstol) == 1 && length(reltol) == 1 - err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate - else # else perform component-wise computations for error - err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate - end + errEstimate = k1 - 2*k2 + k3; + err = getScaledError(y,ynew,errEstimate,h,reltol,abstol,norm) + #if length(abstol) == 1 && length(reltol) == 1 + # err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate + #else + # err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate + #end # check if new solution is acceptable if err <= 1 diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index a43933404..e6ffe2437 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -413,23 +413,9 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, facmin = 1./facmax # maximal step size decrease. ? # in-place calculate xerr./tol + # Get the scaled error and if outside of domain (usually NaN) then make step size smaller by maximum + getScaledError!(x0,xtrial,xerr,reltol,abstol,norm,dof) && return 10., dt*facmin, timout_after_nan - # If reltol and abstol are vectors - # restore old behaviour - if length(abstol) == 1 && length(reltol) == 1 - for d=1:dof - # if outside of domain (usually NaN) then make step size smaller by maximum - isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol + max(abs(x0[d]), abs(xtrial[d]))*reltol) # Eq 4.10 - end - else - # else perform component-wise computations - for d=1:dof - # if outside of domain (usually NaN) then make step size smaller by maximum - isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10 - end - end err = norm(xerr, 2) # Eq. 4.11 newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified if timeout>0 From d68c1a066a0b4d8f378bbac1f5fd5347cfdf9548 Mon Sep 17 00:00:00 2001 From: sonVishal Date: Sat, 9 Apr 2016 13:30:43 +0200 Subject: [PATCH 13/14] Added abs() function which can handle arrays of CompSol to interface-test.jl --- test/interface-tests.jl | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/test/interface-tests.jl b/test/interface-tests.jl index dbb811dcc..816d6773c 100644 --- a/test/interface-tests.jl +++ b/test/interface-tests.jl @@ -45,19 +45,42 @@ ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) ################################################################################ - +# TODO: test a vector-like state variable, i.e. one which can be indexed. +function getZeroVector(y::Array{CompSol,1}) + numElems = length(y) + tmp = Array{CompSol,1}(numElems) + for i = 1:numElems + tmp[i] = CompSol(complex(zeros(2,2)), 0., 0.) + end + return tmp +end +Base.zeros(y::Array{CompSol,1}) = getZeroVector(y) + +function getAbsVector(y::Array{CompSol,1}) + numElems = length(y) + absVec = Array{Float64}(numElems) + for i = 1:numElems + absVec[i] = norm(y[i], 2.) + end + return absVec +end +Base.abs(y::Array{CompSol,1}) = getAbsVector(y) + + +################################################################################ + # define RHSs of differential equations # delta, V and g are parameters function rhs(t, y, delta, V, g) H = [[-delta/2 V]; [V delta/2]] - + rho_dot = -im*H*y.rho + im*y.rho*H x_dot = y.p p_dot = -y.x - + return CompSol( rho_dot, x_dot, p_dot) end - + # inital conditons rho0 = zeros(2,2); rho0[1,1]=1.; @@ -77,6 +100,3 @@ for solver in solvers @test norm(y1[end]-y2[end])<0.1 end println("ok.") - -################################################################################ -# TODO: test a vector-like state variable, i.e. one which can be indexed. From f256c4d980e5dcf793241a7aab949afbefe4971f Mon Sep 17 00:00:00 2001 From: Vishal Sontakke Date: Sun, 18 Sep 2016 17:09:46 +0200 Subject: [PATCH 14/14] Change ROBER testcase relTol and absTol to arrays --- test/runtests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index b66073abc..141390dad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -75,7 +75,7 @@ let ydot end t = [0., 1e11] - t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-9, reltol=1e-9, + t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-9*ones(3), reltol=1e-9*ones(3), maxstep=1e11/10, minstep=1e11/1e18) refsol = [0.2083340149701255e-07, @@ -83,6 +83,7 @@ let 0.9999999791665050] # reference solution at tspan[2] @test norm(refsol-y[end], Inf) < 2e-10 end + include("interface-tests.jl") println("All looks OK")