WIP: Rosenbrock(-W) methods

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.3
 Polynomials
 Compat 0.4.1
+MatrixColorings

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<N
+            ui2 = y[iu+2]
+            vi2 = y[iv+2]
+        else # on last step do BC
+            ui2 = ubc2()
+            vi2 = vbc2()
+        end
+    end
+    return nothing
+end
+fn!(;T_::Type=T, dof_=dof) = zeros(T_,dof_)
+function fn(t, y, ydot)
+    res = fn!()
+    fn!(res, y, ydot)
+    return res
+end
+function fn_ex(t, y)
+    res = fn!()
+    fn!(res, y, 0*y)
+    return res
+end
+
+
+
+# Jacobian is a banded matrix with upper and lower bandwidth 2,
+# c.f. W&H p.148, but Julia does not support banded matrices yet.
+# Thus use a sparse matrix.
+
+function jac!(dfdy, y, ydot, a)
+    # The Jacobian of fn: = df/dy + a* df/dydot
+    ii = 1 # direct index into dfdy.nzval
+    for i=1:N # iterate over double-columns
+        iu, iv = inds(i)
+        ui = y[iu]
+        vi = y[iv]
+        # iu: column belonging to d/du(i)
+        if iu-2>=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

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)