diff --git a/REQUIRE b/REQUIRE index 8e8352dde..2f9450386 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,3 +1,4 @@ julia 0.3 Polynomials Compat 0.4.1 +MatrixColorings diff --git a/examples/bruss1d.jl b/examples/bruss1d.jl new file mode 100644 index 000000000..582005496 --- /dev/null +++ b/examples/bruss1d.jl @@ -0,0 +1,205 @@ +using ODE,DASSL + +const T = Float64 # the datatype used, probably Float64 +Tarr = SparseMatrixCSC +const N = 500 # active gridpoints +const dof = 2N # degrees of freedom (boundary points are not free) +dx = 1/(N+1) +x = dx:dx:1-dx +dae = 0 # index of DAE, ==0 for ODE + +# parameters +alpha = 1/50 +const gamma = alpha/dx^2 + +# BC +const one_ = one(T) +const three = 3*one_ +ubc1() = one_ +ubc2() = one_ +vbc1() = three +vbc2() = three +#IC +u0(x) = 1 + 0.5*sin(2π*x) +v0(x) = 3 * ones(length(x)) + +inds(i) = (2i-1, 2i) +function fn!(res, y, dydt) + # the ode function res = f(y,dydt) + # + # Note the u and v are staggered: + # u = y[1:2:end], v=y[2:2:end] + + # setup for first step + i = 1 + iu, iv = inds(i) + ui0 = ubc1(); vi0 = vbc1() # BC + ui1 = y[iu]; vi1 = y[iv] + ui2 = y[iu+2]; vi2 = y[iv+2] + while i<=N + res[iu] = 1 + ui1^2*vi1 - 4*ui1 + gamma*(ui0 - 2*ui1 + ui2) - dydt[iu] + res[iv] = 3*ui1 - ui1^2*vi1 + gamma*(vi0 - 2*vi1 + vi2) - dydt[iv] + # setup for next step: + i += 1 + iu, iv = inds(i) + ui0 = ui1; vi0 = vi1 + ui1 = ui2; vi1 = vi2 + if i=1 + # du'(i-1)/du(i) + dfdy.nzval[ii] = gamma + ii +=1 + end + if iv-2>=1 + # dv'(i-1)/du(i)==0 + dfdy.nzval[ii] = 0 + ii +=1 + end + # du'(i)/du(i) + dfdy.nzval[ii] = 2*ui*vi - 4 -2*gamma - a + ii += 1 + # dv'(i)/du(i) + dfdy.nzval[ii] = 3 - 2*ui*vi + ii +=1 + if iu+2<=2N + # du'(i+1)/du(i) + dfdy.nzval[ii] = gamma + ii +=1 + end + + # iv: column belonging to d/dv(i) + if iv-2>=1 + # dv'(i-1)/dv(i) + dfdy.nzval[ii] = gamma + ii +=1 + end + # du'(i)/dv(i) + dfdy.nzval[ii] = ui^2 + ii +=1 + # dv'(i)/dv(i) + dfdy.nzval[ii] = -ui^2 -2*gamma - a + ii +=1 + if iu+2<=2N + # du'(i+1)/dv(i)==0 + dfdy.nzval[ii] = 0 + ii +=1 + end + if iv+2<=2N + # dv'(i+1)/dv(i) + dfdy.nzval[ii] = gamma + ii +=1 + end + end + return nothing +end +function jac!(;T_::Type=T, dof_=dof) + B = (ones(T_,dof_-2), ones(T_,dof_-1), ones(T_,dof_), ones(T_,dof_-1), ones(T_,dof_-2)) + spdiagm(B, (-2,-1,0,1,2)) +end +function jac(t, y, ydot, a) + J = jac!() + jac!(J, y, ydot, a) + return J +end +function jac_ex(t, y) + return jac(t, y, y, 0) +end + +refsol = T[0.9949197002317599, 3.0213845767604077, 0.9594350193986054, 3.0585989778165419, 0.9243010095428502, 3.0952478919989637, 0.8897959106772672, + 3.1310118289054731, 0.8561653620284367, 3.1656101198770159, 0.8236197147449046, 3.1988043370624344, 0.7923328094811884, 3.2303999530641514, + 0.7624421042573115, 3.2602463873623941, 0.7340499750795348, 3.2882356529108807, 0.7072259700779899, 3.3142998590079271, 0.6820097782458483, + 3.3384078449410937, 0.6584146743834650, 3.3605612157873943, 0.6364312187752559, 3.3807900316323134, 0.6160310186921587, 3.3991483695914764, + 0.5971703941198909, 3.4157099395342736, 0.5797938277687891, 3.4305638938070224, 0.5638371159206763, 3.4438109320334580, 0.5492301695479158, + 3.4555597666485198, 0.5358994429426996, 3.4659239846027008, 0.5237699892215797, 3.4750193162238476, 0.5127671585747183, 3.4829613034792271, + 0.5028179665048467, 3.4898633463634923, 0.4938521662914935, 3.4958350971335204, 0.4858030633656755, 3.5009811668111510, 0.4786081100251151, + 3.5054001059792705, 0.4722093177200750, 3.5091836216744015, 0.4665535216425440, 3.5124159935026285, 0.4615925290790646, 3.5151736544621075, + 0.4572831793403656, 3.5175249049438184, 0.4535873393501199, 3.5195297317024448, 0.4504718553589467, 3.5212397070273984, 0.4479084778719241, + 3.5226979467564341, 0.4458737738041973, 3.5239391090719634, 0.4443490371324889, 3.5249894191569453, 0.4433202068820853, 3.5258667077466495, + 0.4427777991494095, 3.5265804544017270, 0.4427168579654424, 3.5271318289682063, 0.4431369281018266, 3.5275137272135266, 0.4440420513508381, + 3.5277107990730161, 0.4454407863109616, 3.5276994703501980, 0.4473462502188303, 3.5274479611304068, 0.4497761798232572, 3.5269163066324394, + 0.4527530066369863, 3.5260563887768472, 0.4563039400688689, 3.5248119894251024, 0.4604610498812091, 3.5231188790654930, 0.4652613370907894, + 3.5209049576992761, 0.4707467798082714, 3.5180904678044698, 0.4769643375804777, 3.5145883024867057, 0.4839658945842979, 3.5103044351908528, + 0.4918081185812277, 3.5051385005173827, 0.5005522089940899, 3.4989845585737802, 0.5102635039989190, 3.4917320776245013, 0.5210109134090777, + 3.4832671712209993, 0.5328661417420937, 3.4734741260299615, 0.5459026646938675, 3.4622372546582585, 0.5601944229089820, 3.4494431032230182, + 0.5758142001453760, 3.4349830354873361, 0.5928316594749734, 3.4187562033108012, 0.6113110218368440, 3.4006728962523969, 0.6313083867734524, + 3.3806582409098729, 0.6528687160104193, 3.3586561928427350, 0.6760225267555723, 3.3346337311179157, 0.7007823726569721, 3.3085851288057930, + 0.7271392249346637, 3.2805361342349380, 0.7550589020044152, 3.2505478606008622, 0.7844787296769868, 3.2187201496972175, 0.8153046416214843, + 3.1851941538893653, 0.8474089465959840, 3.1501538739882800, 0.8806289904192589, 3.1138264039027113, 0.9147669230929857, 3.0764806689389470, + 0.9495907429372025, 3.0384245041548366, 0.9848367306701233] # reference sol from Hairer has 1.36 scd +refsolinds = collect(1:7:dof) # refsol is only at components 1:7:2N + + +ic = zeros(T,2N) +ic[1:2:2N] = u0(x) +ic[2:2:2N] = v0(x) +tspan = T[0.0, 10.0] # integration interval +tspan0 = T[0.0, 0.001] # integration interval +tspan1 = T[0.0, 0.3] # integration interval + +#t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac) + +## ode23s +tout, yout = ode23s(fn_ex, ic, tspan0, jacobian=jac_ex) +@time tout, yout = ode23s(fn_ex, ic, tspan, jacobian=jac_ex) +@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf) + +# ## DASSL +t,ydassl = dasslSolve(fn, ic, tspan0, jacobian=jac) +@time t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac) +@show norm(abs(ydassl[end][refsolinds]-refsol)./refsol, Inf) + +## ROSW +# for some reason the abstol and reltol need to be way lower to reach +# the same precision. Presumably an error in the step control. +abstol, reltol = 1e-2, 1e-3 +# No jac +tout, yout = ode_rosw(fn!, ic, tspan0) +@time tout, yout = ode_rosw(fn!, ic, tspan, abstol=abstol, reltol=reltol); +@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf) + +# Analytic Jacobian, same as previous +tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=(jac!, jac!())) +@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=(jac!, jac!()), abstol=abstol, reltol=reltol); +@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf) + +# Numerical colored Jacobian +tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=jac!()) +@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=jac!(), abstol=abstol, reltol=reltol); +@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf) + +nothing diff --git a/examples/hb1dae.jl b/examples/hb1dae.jl new file mode 100644 index 000000000..e6902685f --- /dev/null +++ b/examples/hb1dae.jl @@ -0,0 +1,87 @@ +# hb1dae reformulates the hb1ode example as a differential-algebraic equation (DAE) problem. +# See Matlab help and http://radio.feld.cvut.cz/matlab/techdoc/math_anal/ch_8_od8.html#670396 + +using ODE, DASSL +M = diagm([1.,1,1]) +M[3,3] = 0. + +function hb1dae(t, y,ydot) + res = zeros(length(y)) + hb1dae!(res, y,ydot) + return res +end + +function hb1dae!(res, y,ydot) + res[1] = -0.04*y[1] + 1e4*y[2]*y[3] - ydot[1] + res[2] = 0.04*y[1] - 1e4*y[2]*y[3] - 3e7*y[2]^2 - ydot[2] + res[3] = y[1] + y[2] + y[3] - 1 + return nothing +end +function Jhb1dae!(res, y, ydot, a) + res[1,1] = -0.04 - a + res[1,2] = 1e4*y[3] + res[1,3] = 1e4*y[2] + + res[2,1] = 0.04 + res[2,2] = -1e4*y[3] - 3e7*2*y[2] -a; + res[2,3] = -1e4*y[2] + + res[3,1] = 1. + res[3,2] = 1. + res[3,3] = 1. + + return nothing +end +function Jhb1dae(t, y, ydot, a) + jac = zeros(length(y),length(y)) + Jhb1dae!(jac, y, ydot, a) + return jac +end + + +tspan = [0, 4e6] +y0 = [1., 0, 0] + +t,ydassl = dasslSolve(hb1dae, y0, tspan, jacobian=Jhb1dae) +@time t,ydassl = dasslSolve(hb1dae, y0, tspan, jacobian=Jhb1dae) +y = hcat(ydassl...); + +## ROSW +t,yrosw = ode_rosw(hb1dae!, y0, tspan) #, jacobian=Jhb1dae!) +@time t,yrosw = ode_rosw(hb1dae!, y0, tspan)#, jacobian=Jhb1dae!) +yr = hcat(yrosw...); + +println("Relative difference between DASSL vs ROSW, no Jac:") +println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end])) + +# println("Using fixed step ROSW") +# tspan = t +# t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan) +# @time t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan) +# println("Relative difference between DASSL vs fixed step, no Jac:") +# println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end])) + +# with Jacobian +tspan = [0, 4e6] +t,yrosw = ode_rosw(hb1dae!, y0, tspan, jacobian=Jhb1dae!) +@time t,yrosw = ode_rosw(hb1dae!, y0, tspan, jacobian=Jhb1dae!) +yr = hcat(yrosw...); + +println("Relative difference between DASSL vs ROSW with Jac:") +println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end])) + + +println("Using ros_rodas3") +t,yrosw = ODE.oderosw_adapt(hb1dae!, y0, tspan, ODE.bt_ros_rodas3) +@time t,yrosw = ODE.oderosw_adapt(hb1dae!, y0, tspan, ODE.bt_ros_rodas3) +println("Relative difference between DASSL vs Robas2, no Jac:") +println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end])) + + + + +# using Winston +# plot(t, y[1,:], xlog=true) +# oplot(t, y[2,:]*1e4, xlog=true) +# oplot(t, y[3,:], xlog=true) +# oplot(t, y[1,:], xlog=true) diff --git a/examples/van_der_pol_rosw.jl b/examples/van_der_pol_rosw.jl new file mode 100644 index 000000000..2bdb32c28 --- /dev/null +++ b/examples/van_der_pol_rosw.jl @@ -0,0 +1,193 @@ +# Van der Pol using DASSL.jl +using ODE +using DASSL +abstol = 1e-6 +reltol = 1e-6 + +const mu = 100. + +# Explicit equation +function vdp(t, y) + # vdp(t, y) + # + # ydot of van der Pol equation. (consider using rescaled eq. so + # period stays the same for different mu. Cf. Wanner & Hairer 1991, + # p.5) + return [y[2], mu^2*((1-y[1]^2) * y[2] - y[1])] +end +function Jvdp(t, y) + # Jvdp(t, y, mu) + # + # Jacobian of vdp + return [0 1; + -mu^2*(y[2]*2*y[1] + 1) mu^2*(1-y[1]^2)] +end + +# implicit equation +function vdp_impl(y, ydot) + # vdp_impl(t, y) + # + # Implicit van der Pol equation. (consider using rescaled eq. so + # period stays the same for different mu. Cf. Wanner & Hairer 1991, + # p.5) + return vdp(0,y) - ydot +end +function Jvdp_impl(y, ydot, α) + # d vdp_impl /dy + α d vdp_impl/dydot + # + # Jacobian + return Jvdp(0,y) - α * eye(2) +end + +# Results 24dfdc9cacfd6 +# elapsed time: 0.025767095 seconds (25194952 bytes allocated) +# elapsed time: 0.00997568 seconds (6674864 bytes allocated) +# Adaptive step: abs error of ode23s vs ref: +# [0.012513592025336528,0.013225223452835055] +# Adaptive step: abs error of ode_rosw vs ref: +# [0.012416936355840624,0.013124961519338618] + + +# implicit inplace equation +function vdp_impl!(res, y, ydot) + # vdp_impl(t, y) + # + # Implicit van der Pol equation. (consider using rescaled eq. so + # period stays the same for different mu. Cf. Wanner & Hairer 1991, + # p.5) + # + # Note, that ydot can be used as res here. + res[1] = y[2] - ydot[1] + res[2] = mu^2*((1-y[1]^2) * y[2] - y[1]) - ydot[2] + return nothing +end +function Jvdp_impl!(res, y, ydot, α) + # d vdp_impl /dy + α d vdp/dydot + # + # Jacobian + + res[1,1] = 0 - α + res[1,2] = 1 + res[2,1] = -mu^2*(y[2]*2*y[1] + 1) + res[2,2] = mu^2*(1-y[1]^2) - α + return nothing +end + +# elapsed time: 0.02611999 seconds (25194952 bytes allocated) +# elapsed time: 0.006630583 seconds (4025352 bytes allocated) +# elapsed time: 0.005866466 seconds (2670408 bytes allocated) (newer version) +# Adaptive step: abs error of ode23s vs ref: +# [0.012513592025336528,0.013225223452835055] +# Adaptive step: abs error of ode_rosw vs ref: +# [0.012416936355840624,0.013124961519338618] + + +### +# reference as t=2 +## +refsol = [0.1706167732170483e1, -0.8928097010247975e0] +y0 = [2.,0.]; +tstart = 0. +tend = 2. + + + +### +# Fixed step +### +# nt = 50_000; +# tspan = linspace(tstart, tend, nt) + +# t,yout1 = ode4s(vdp, y0, tspan; jacobian=Jvdp) +# @time t,yout1 = ode4s(vdp, y0, tspan; jacobian=Jvdp) + +# t,yout2 = ode_rosw_fixed(vdp_impl, Jvdp_impl, y0, tspan) +# @time t,yout2 = ode_rosw_fixed(vdp_impl, Jvdp_impl, y0, tspan) + +# println("Fixed step: abs error of ode4s vs ref:") +# println(yout1[end]-refsol) +# println("Fixed step: abs error of ode_rosw vs ref:") +# println(yout2[end]-refsol) + +# #### adaptive +tspan = linspace(tstart, tend, 2) + +t,yout3 = ode23s(vdp, y0, [0,0.1]; jacobian=Jvdp) +gc() +@time t,yout3 = ode23s(vdp, y0, tspan; jacobian=Jvdp) + +t,yout4 = ode_rosw(vdp_impl!, y0, [0, 0.1], jacobian=Jvdp_impl!) +gc() +@time t,yout4 = ode_rosw(vdp_impl!, y0, tspan, jacobian=Jvdp_impl!) + +# numerical jacobian +t,yout5 = ode_rosw(vdp_impl!, y0, [0, 0.1]) +gc() +@time t,yout5 = ode_rosw(vdp_impl!, y0, tspan) + + +println("Adaptive step: abs error of ode23s vs ref:") +println(yout3[end]-refsol) +println("Adaptive step: abs error of ode_rosw vs ref:") +println(yout4[end]-refsol) +println("Adaptive step: ROSW abs error with numerical Jacobian:") +println(yout5[end]-refsol) + + +# #################### +# # DAE: reduced Van der Pol + + + +# # implicit inplace equation +# function dae_vdp_impl!(res, y, ydot) +# # dae_vdp_impl(t, y) +# # +# # Implicit, reduced van der Pol equation. mu->\infinity + +# res[1] = y[2] - ydot[1] +# res[2] = y[1] - (y[2]^3/3 - y[2]) +# return nothing +# end +# function Jdae_vdp_impl!(res, y, ydot, α) +# # d dae_vdp_impl /dy + α d dae_vdp_impl/dydot +# # +# # Jacobian + +# res[1,1] = 0 - α +# res[1,2] = 1 +# res[2,1] = 1 +# res[2,2] = -(y[2]^2 - 1) +# return nothing +# end + +# t,yout5 = ode_rosw(dae_vdp_impl!, Jdae_vdp_impl!, y0, [0, 0.1]) +# gc() +# @time t,yout5 = ode_rosw(dae_vdp_impl!, Jdae_vdp_impl!, y0, tspan) + + + +# # implicit inplace equation +# function dae_vdp(t, y, ydot) +# # dae_vdp_impl(t, y) +# # +# # Implicit, reduced van der Pol equation. mu->\infinity + +# [y[2] - ydot[1], y[1] - (y[2]^3/3 - y[2])] +# end +# function Jdae_vdp(t, y, ydot, α) +# # d dae_vdp_impl /dy + α d dae_vdp_impl/dydot +# # +# # Jacobian +# res = zeros(2,2) +# res[1,1] = 0 - α +# res[1,2] = 1 +# res[2,1] = 1 +# res[2,2] = -(y[2]^2 - 1) +# return res +# end + +# t,ydassl = dasslSolve(dae_vdp, y0, tspan) + + + diff --git a/src/ODE.jl b/src/ODE.jl index a364988d8..6738f12ed 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,6 +62,7 @@ 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) + # 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) @@ -84,7 +85,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) pow = -(2. + log10(max(d1, d2)))/(p + 1.) h1 = 10.^pow end - return tdir*min(100*h0, h1), tdir, f0 + return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0 end # isoutofdomain takes the state and returns true if state is outside @@ -423,4 +424,7 @@ const ms_coefficients4 = [ 1 0 0 0 -9/24 37/24 -59/24 55/24] +## Rosenbrock-Wanner methods +include("rosw.jl") + end # module ODE diff --git a/src/rosw.jl b/src/rosw.jl new file mode 100644 index 000000000..fa4bd1b7b --- /dev/null +++ b/src/rosw.jl @@ -0,0 +1,509 @@ +# Rosenbrock-Wanner methods +########################### +# +# Main references: +# - Wanner & Hairer 1996 +# - PETSc http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSROSW.html + +using MatrixColorings +export ode_rosw, ode_rosw_fixed, ode_rodas3 + +# TODO: +# - Jacobian coloring: https://github.com/mlubin/ReverseDiffSparse.jl +# http://wiki.cs.purdue.edu:9835/coloringpage/abstracts/acyclic-SISC.pdf +# - AD Jacobian: https://github.com/mlubin/ReverseDiffSparse.jl +# - fix fixed step solver + +# Rosenbrock-W methods are typically specified for autonomous DAE: +# Mẋ = f(x) +# +# by the stage equations +# M kᵢ = hf(x₀ + + Σⱼ aᵢⱼkⱼ) + h J Σⱼ γᵢⱼkⱼ +# +# and step completion formula +# x₁ = x₀ + \Sigmaⱼ bⱼ kⱼ +# +# The method used here uses transformed equations as done in PETSc. + + +# rosw_poststep(xtrial, err, t, dt, steps) = (showcompact(t), print(", ") , +# showcompact(dt), print(", "), +# showcompact(err), print(", "), +# println(steps)) + +rosw_poststep(xtrial, err, t, dt, steps) = nothing + +# Tableaus +########## + +immutable TableauRosW{Name, S, T} <: Tableau{Name, S, T} + order::(@compat(Tuple{Vararg{Int}})) # the order of the methods + a::Matrix{T} + γ::Matrix{T} + # one or several row vectors. First row is used for the step, + # second for error calc. + b::Matrix{T} + function TableauRosW(order,a,γ,b) + @assert isa(S,Integer) + @assert isa(Name,Symbol) + @assert S==size(γ,1)==size(a,2)==size(γ,1)==size(a,2)==size(b,2) + @assert size(b,1)==length(order) + new(order,a,γ,b) + end +end +function TableauRosW{T}(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), + a::Matrix{T}, γ::Matrix{T}, b::Matrix{T}) + TableauRosW{name,size(b,2),T}(order, a, γ, b) +end +function TableauRosW(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), T::Type, + a::Matrix, γ::Matrix, b::Matrix) + TableauRosW{name,size(b,2),T}(order, convert(Matrix{T},a), + convert(Matrix{T},γ), convert(Matrix{T},b) ) +end +conv_field{T,N}(D,a::Array{T,N}) = convert(Array{D,N}, a) +function Base.convert{Tnew<:Real,Name,S,T}(::Type{Tnew}, tab::TableauRosW{Name,S,T}) + # Converts the tableau coefficients to the new type Tnew + newflds = () + @compat for n in fieldnames(tab) + fld = getfield(tab,n) + if eltype(fld)==T + newflds = tuple(newflds..., conv_field(Tnew, fld)) + else + newflds = tuple(newflds..., fld) + end + end + TableauRosW{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way? +end + + + +# Transformed Tableau, used only internally +immutable TableauRosW_T{Name, S, T} <: Tableau{Name, S, T} + order::(@compat(Tuple{Vararg{Int}})) # the order of the methods + a::Matrix{T} # this is TableauRosW.a transformed + γinv::Matrix{T} + b::Matrix{T} # this is TableauRosW.b transformed + # derived quantities: + γii::T + c::Matrix{T} # = (tril(btab.γinv)-diagm(diag(btab.γinv))) +end +function tabletransform{Name,S,T}(rt::TableauRosW{Name,S,T}) + # the code only works if + if !all(x->x==rt.γ[1],diag(rt.γ)) + error("This Rosenbrock implementation only works for tableaus with γ_ii==γ_jj for all i,j.") + end + γii = rt.γ[1,1] + γinv = inv(rt.γ) + ahat = rt.a * γinv + bhat = similar(rt.b) + bhat[1,:] = squeeze(rt.b[1,:]*γinv,1) + bhat[2,:] = squeeze(rt.b[2,:]*γinv,1) + c = (tril(γinv)-diagm(diag(γinv))) # negative of W&H definition + return TableauRosW_T{Name,S,T}(rt.order, ahat, γinv, bhat, γii, c) +end + + +## tableau for ros34pw2 http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSROSWRA34PW2.html +const bt_ros34pw2 = TableauRosW(:ros34pw2, (3,4), Float64, + [0 0 0 0 + 8.7173304301691801e-01 0 0 0 + 8.4457060015369423e-01 -1.1299064236484185e-01 0 0 + 0 0 1. 0], + [4.3586652150845900e-01 0 0 0 + -8.7173304301691801e-01 4.3586652150845900e-01 0 0 + -9.0338057013044082e-01 5.4180672388095326e-02 4.3586652150845900e-01 0 + 2.4212380706095346e-01 -1.2232505839045147e+00 5.4526025533510214e-01 4.3586652150845900e-01], + [2.4212380706095346e-01 -1.2232505839045147e+00 1.5452602553351020e+00 4.3586652150845900e-01 + 3.7810903145819369e-01 -9.6042292212423178e-02 5.0000000000000000e-01 2.1793326075422950e-01] + ) + +## tableau for RODAS3 http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSROSWRODAS3.html +const bt_ros_rodas3 = TableauRosW(:ros_rodas3, (3,4), Rational{Int}, + [0 0 0 0 + 0 0 0 0 + 1 0 0 0 + 3//4 -1//4 1//2 0], + [1//2 0 0 0 + 1 1//2 0 0 + -1//4 -1//4 1//2 0 + 1//12 1//12 -2//3 1//2], + [5//6 -1//6 -1//6 1//2 + 3//4 -1//4 1//2 0] + ) + +################### +# Fixed step solver +################### +ode_rosw_fixed(fn!, x0, tspan; kwargs...) = oderosw_fixed(fn!, x0, tspan, bt_ros34pw2, kwargs...) +function oderosw_fixed{N,S}(fn!, x0::AbstractVector, tspan, + btab::TableauRosW{N,S}; + jacobian=numerical_jacobian(fn!, 1e-6, 1e-6) + ) + # TODO: refactor with oderk_fixed + Et, Exf, Tx, btab = make_consistent_types(fn!, x0, tspan, btab) + btab = tabletransform(btab) + dof = length(x0) + + xs = Array(Tx, length(tspan)) + allocate!(xs, x0, dof) + xs[1] = deepcopy(x0) + + tspan = convert(Vector{Et}, tspan) + # work arrays: + k = Array(Tx, S) # stage variables + allocate!(k, x0, dof) + ks = zeros(Exf,dof) # work vector for one k + jac_store = zeros(Exf,dof,dof) # Jacobian storage + u = zeros(Exf,dof) # work vector + udot = zeros(Exf,dof) # work vector + + # allocate!(ks, x0, dof) # no need to allocate as fn! is not in-place + xtmp = similar(x0, Exf, dof) + for i=1:length(tspan)-1 + dt = tspan[i+1]-tspan[i] + rosw_step!(xs[i+1], fn!, jacobian, xs[i], dt, dof, btab, + k, jac_store, ks, u, udot) + end + return tspan, xs +end + + +ode_rosw(fn!, x0, tspan;kwargs...) = oderosw_adapt(fn!, x0, tspan, bt_ros34pw2; kwargs...) +ode_rodas3(fn!, x0, tspan;kwargs...) = oderosw_adapt(fn!, x0, tspan, bt_ros_rodas3; kwargs...) +function oderosw_adapt{N,S}(fn!, x0::AbstractVector, tspan, btab::TableauRosW{N,S}; + reltol = 1.0e-5, abstol = 1.0e-8, + norm=Base.norm, + minstep=abs(tspan[end] - tspan[1])/1e18, + maxstep=abs(tspan[end] - tspan[1])/2.5, + initstep=0., + jacobian=numerical_jacobian(fn!, reltol, abstol) +# points=:all + ) + # TODO: refactor with oderk_adapt + + # FIXME: add interpolation + if length(tspan)>2 + warn("specified output times not supported yet") + tspan = [tspan[1], tspan[end]] + points=:all + else + points=:all + end + ## Figure types + fn_expl = (t,x)->(out=similar(x); fn!(out, x, x*0); out) + Et, Exf, Tx, btab = make_consistent_types(fn!, x0, tspan, btab) + btab = tabletransform(btab) + # parameters + order = minimum(btab.order) + timeout_const = 5 # after step reduction do not increase step for + # timeout_const steps + + ## Setup + dof = length(x0) + tspan = convert(Vector{Et}, tspan) + t = tspan[1] + tstart = tspan[1] + tend = tspan[end] + + # work arrays: + x = similar(x0, Exf, dof) # x at time t (time at beginning of step) + x[:] = x0 # fill with IC + xtrial = similar(x0, Exf, dof) # trial solution at time t+dt + xerr = similar(x0, Exf, dof) # error of trial solution + k = Array(Tx, S) # stage variables + allocate!(k, x0, dof) + ks = zeros(Exf,dof) # work vector for one k + if isa(jacobian, Tuple) # Jacobian function and sparse pattern/storage provided + jac_store = jacobian[2] + jacobian = jacobian[1] + elseif isa(jacobian, SparseMatrixCSC) # Sparse storage provided, make colored Jacobian (TODO) + jac_store = jacobian + jacobian = numerical_jacobian(fn!, reltol, abstol, jac_store) + else # nothing provided, make dense numerical jacobian + jac_store = zeros(Exf,dof,dof) # Jacobian storage + end + u = zeros(Exf,dof) # work vector + udot = zeros(Exf,dof) # work vector + + # output xs + nsteps_fixed = length(tspan) # these are always output + xs = Array(Tx, nsteps_fixed) + allocate!(xs, x0, dof) + xs[1] = x0 + + # Option points determines where solution is returned: + if points==:all + tspan_fixed = tspan + tspan = Et[tstart] + iter_fixed = 2 # index into tspan_fixed + sizehint!(tspan, nsteps_fixed) + elseif points!=:specified + error("Unrecognized option points==$points") + end + + # Time + dt, tdir, f0 = hinit(fn_expl, x, tstart, tend, order, reltol, abstol) + if initstep!=0 + dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") + end + # Diagnostics + dts = Et[] + errs = Float64[] + steps = [0,0] # [accepted, rejected] + + ## Integration loop + islaststep = abs(t+dt-tend)<=eps(tend) ? true : false + timeout = 0 # for step-control + iter = 2 # the index into tspan and xs + while true + rosw_step!(xtrial, fn!, jacobian, x, dt, dof, btab, + k, jac_store, ks, u, udot) + # Completion again for embedded method, see line 927 of + # http://www.mcs.anl.gov/petsc/petsc-current/src/ts/impls/rosw/rosw.c.html#TSROSW + xerr[:] = 0.0 + for s=1:S + for d=1:dof + xerr[d] += (btab.b[2,s]-btab.b[1,s])*k[s][d] + end + end + err, newdt, timeout = stepsize_hw92!(dt, tdir, x, xtrial, xerr, order, timeout, + dof, abstol, reltol, maxstep, norm) + + rosw_poststep(xtrial, err, t, dt, steps) + if err<=1.0 # accept step + # diagnostics + steps[1] +=1 + push!(dts, dt) + push!(errs, err) + + # Output: + # f0 = k[1] + # f1 = fn_expl(t+dt, xtrial) + if points==:specified + # interpolate onto given output points + while iter-1= tdir*tend + dt = tend-t + islaststep = true # next step is the last, if it succeeds + end + elseif abs(newdt)= tdir*tend dt = tend-t - laststep = true # next step is the last, if it succeeds + islaststep = true # next step is the last, if it succeeds end elseif abs(newdt)0 no step size increase is allowed, timeout is # decremented in here. # - # Returns the error, newdt and the number of timeout-steps + # Returns the error, newdt and the number of timeout-steps. + # Updates xerr in-place with the relative error. # # TODO: # - allow component-wise reltol and abstol? @@ -464,7 +465,6 @@ function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1) # # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 ) # this is O(3). TODO for higher order. - theta = (tquery-t)/dt for i=1:length(y0) y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) * diff --git a/test/implicit-eqs.jl b/test/implicit-eqs.jl new file mode 100644 index 000000000..ada0eb754 --- /dev/null +++ b/test/implicit-eqs.jl @@ -0,0 +1,78 @@ +using ODE +using Base.Test + +tol = 1e-2 + +solvers = [ + ## Stiff + ODE.ode_rosw + ] + +for solver in solvers + println("using $solver") + # dy + # -- = 6 ==> y = 6t + # dt + f = (res,y,ydot)->(res[:]=6.0-ydot; nothing) + t,y=solver(f, [0.], [0.,1]) + y = vcat(y...) + @test maximum(abs(y-6t)) < tol + + # dy + # -- = 2t ==> y = t.^2 + # dt + # t,y=solver((t,y)->2t, 0., [0:.001:1;]) + # @test maximum(abs(y-t.^2)) < tol + + + # dy + # -- = y ==> y = y0*e.^t + # dt + f = (res,y,ydot)->(res[:]=y-ydot;nothing) + t,y=solver(f, [1.], [0,1]) + y = vcat(y...) + @test maximum(abs(y-e.^t)) < tol + + t,y=solver((res,y,ydot)->(res[:]=y-ydot), [1.], [1, 0]) + y = vcat(y...) + @test maximum(abs(y-e.^(t-1))) < tol + + # dv dw + # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t) + # dt dt + # + # y = [v, w] + + function ff(res,y,ydot) + res[1] = -y[2]-ydot[1] + res[2] = y[1] -ydot[2] + end + + t,y=solver(ff, [1., 2.], [0, 2*pi]) + ys = hcat(y...).' # convert Vector{Vector{Float}} to Matrix{Float} + @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol +end + +# rober testcase from http://www.unige.ch/~hairer/testset/testset.html +let + println("ROBER test case") + function f(t, y) + ydot = similar(y) + ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3] + ydot[3] = 3.0e7*y[2]*y[2] + ydot[2] = -ydot[1] - ydot[3] + ydot + end + f!(res, y, ydot) = (res[:] = f(0, y) - ydot; nothing) + + t = [0., 1e11] + t,y = ode_rosw(f!, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, + 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) < 3e-10 +end + +println("All looks OK") diff --git a/test/runtests.jl b/test/runtests.jl index 64f57e0b7..9376f8d8d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,6 +23,7 @@ solvers = [ # fixed-step ODE.ode4s_s, ODE.ode4s_kr, +# ODE.ode_rosw, # adaptive ODE.ode23s] @@ -80,5 +81,6 @@ let @test norm(refsol-y[end], Inf) < 2e-10 end include("interface-tests.jl") +include("implicit-eqs.jl") println("All looks OK")