basic implementation working

Author: vyudu

lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl

@@ -7,7 +7,7 @@ using SymbolicIndexingInterface
 using LinearAlgebra
 using ADTypes
 using UnPack
-import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, initialize!, perform_step!
+import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required
 import CommonSolve
 import DifferentiationInterface as DI
 
@@ -25,6 +25,7 @@ end
 
 include("cache.jl")
 include("solve.jl")
+include("alg_utils.jl")
 
 export IDSolve
 

lib/ImplicitDiscreteSolve/src/alg_utils.jl

@@ -0,0 +1,19 @@
+function SciMLBase.isautodifferentiable(alg::IDSolve)
+    true
+end
+function SciMLBase.allows_arbitrary_number_types(alg::IDSolve)
+    true
+end
+function SciMLBase.allowscomplex(alg::IDSolve)
+    true
+end
+
+SciMLBase.isdiscrete(alg::IDSolve) = true
+
+isfsal(alg::IDSolve) = false
+alg_order(alg::IDSolve) = 0
+beta2_default(alg::IDSolve) = 0
+beta1_default(alg::IDSolve, beta2) = 0
+
+dt_required(alg::IDSolve) = false
+isdiscretealg(alg::IDSolve) = true

lib/ImplicitDiscreteSolve/src/solve.jl

@@ -5,7 +5,6 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     @unpack state, prob = cache
     state.u .= uprev
     state.t_next = t
-    @show state
     prob = remake(prob, p = state)
 
     u = solve(prob, nlsolve)
@@ -32,15 +31,3 @@ function initialize!(integrator, cache::IDSolveCache)
     end
     cache.prob = prob
 end
-
-#### unnecessary
-#function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg) 
-#    f = prob.f
-#    t_i = prob.tspan[1]
-#    u0 = state_values(prob)
-#    p = parameter_values(prob)
-#
-#    _f(resid, u_next, p) = f(resid, u_next, p.u, p.p, p.t)
-#    state = ImplicitDiscreteState(u0, p, t_i)
-#    nlprob = NonlinearProblem(_f, u0, state)
-#end

lib/ImplicitDiscreteSolve/test/runtests.jl

@@ -6,21 +6,24 @@ using SimpleNonlinearSolve
 
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Discrete System" begin
-    function lotkavolterra(u, p, t) 
+    function lotkavolterra(u, p, t)
         [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]]
     end
 
     function f!(resid, u_next, u, p, t)
-        @. resid = u_next - 0.01*lotkavolterra(u_next, p, t)
+        lv = lotkavolterra(u_next, p, t)
+        resid[1] = u_next[1] - u[1] - 0.01*lv[1]
+        resid[2] = u_next[2] - u[2] - 0.01*lv[2]
+        nothing
     end
     u0 = [1., 1.]
-    tspan = (0., 1.)
+    tspan = (0., 0.5)
 
     idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01)
     idsol = solve(idprob, IDSolve(SimpleNewtonRaphson()))
 
     oprob = ODEProblem(lotkavolterra, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
 
-    @test idsol[end] ≈ osol[end]
+    @test isapprox(idsol[end], osol[end], atol = 0.01)
 end