Testing lazy W operator

testW.jl

using LinearAlgebra, Random, OrdinaryDiffEq, DiffEqOperators
println("-------------------------")
println("Linear univariate problem")
α = -0.5; u0 = 1.0; tspan = (0.0,1.0)
fun = ODEFunction((u,p,t) -> α*u;
                  analytic=(u0,p,t) -> exp(α*t)*u0)
prob = ODEProblem(fun,u0,tspan)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Finite difference error = $err")
fun = ODEFunction((u,p,t) -> α*u;
                  jac_prototype=DiffEqScalar(α),
                  analytic=(u0,p,t) -> exp(α*t)*u0)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Exact Jacobian error = $err")

println("-------------------------")
println("Linear bivariate problem")
A = [-1.0 0.0; 0.0 -0.5]; u0 = [1.0, 1.0]; tspan = (0.0,1.0)
fun = ODEFunction((u,p,t) -> A*u;
                  analytic=(u0,p,t) -> exp(A*t)*u0)
prob = ODEProblem(fun,u0,tspan)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Finite difference error = $err")
fun = ODEFunction((u,p,t) -> A*u;
                  jac_prototype=DiffEqArrayOperator(A),
                  analytic=(u0,p,t) -> exp(A*t)*u0)
prob = ODEProblem(fun,u0,tspan)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Exact Jacobian error = $err")

println("-------------------------")
println("Linear bivariate problem with mass matrix")
A = [-1.0 0.0; 0.0 -0.5]; u0 = [1.0, 1.0]; tspan = (0.0,1.0)
srand(0); mm = [2.0 0.0; 0.0 1.0]; mmA = [-0.5 0.0; 0.0 -0.5]
fun = ODEFunction((u,p,t) -> A*u;
                  analytic=(u0,p,t) -> exp(mmA*t)*u0)
prob = ODEProblem(fun,u0,tspan; mass_matrix=mm)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Finite difference error = $err")
fun = ODEFunction((u,p,t) -> A*u;
                  jac_prototype=DiffEqArrayOperator(A),
                  analytic=(u0,p,t) -> exp(mmA*t)*u0)
prob = ODEProblem(fun,u0,tspan; mass_matrix=mm)
sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
println("Exact Jacobian error = $err")