Explicit Jacobian on DAEs?


#21

thanks


#22

Is this how you woudl go about it?


using ForwardDiff

function test(u)
  solve(PKPD3_Prob,IDA())[1]#,save_everystep=false)[1]
end
ForwardDiff.jacobian(test,IC)

This jacobiam gives a matrix of zeros, which I am certain is not correct so I know I am doing something wrong, but can’t see what- except that all I am returning is aconstant function with the value of PKPD3 at time zero.

I just wanted to compare the values at time 0 from my anlytical and the ForwardDiff at time 0, but I am clearly doing something wrong


#23

Look at http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1


#24

Thanks again


#25

So been looking at the sensitivity analysis, and I don’t get very well how the jacobian is supposed to work - so I tested some of you code in the page tp see if that would help my understanding and I found that the results between ForwardDiff and DiffResults don’t match -
I just copied and pasted, except for one little typo (I believe) correction noted in the code comment, (before my correction I run your code exactly as it was - and the result stayed the same after my correction):

The Forwarddiff code


function f(du,u,p,t)
  dx = p[1]*u[1] - p[2]*u[1]*u[2]
  dy = -p[3]*u[2] + u[1]*u[2]
end

p = [1.5,1.0,3.0,1.0]
u0 = [1.0;1.0]
prob = ODEProblem(f,u0,(0.0,10.0),p)

function test_f(p)
  _prob = remake(prob;u0=convert.(eltype(p),prob.u0),p=p)
  #solve(prob,Vern9(),save_everystep=false)[end]           # I run this once and corrected below once
  solve(_prob,Vern9(),save_everystep=false)[end]
end

using ForwardDiff
fd_res = ForwardDiff.jacobian(test_f,p)

returns a jacobian of 2x4 of only zeroes, no matter if I run the corrected or uncorrected code

I think this i wrong as


using DiffResults
res = DiffResults.JacobianResult(u0,p) # Build the results object
DiffResults.jacobian!(res,u0) # Populate it with the results
val = DiffResults.value(res) # This is the sol[end]
jac = DiffResults.jacobian(res) # This is dsol/dp

returns

Jacobian = [1 0 0 0; 1 0 0 0]

ForwardDiff is working well as when I try something simpler

function t(x)
    [x[1]+x[2]
    2*x[1]*x[2]]
end

fd_res = ForwardDiff.jacobian(t,[1.0,1])

it returns teh correct jacobian [1 1;2 2]

It seems to me that test_f§ is sort of returning a constant array instead of a function so all derivatives are 0, but I am not sure if that assesment it right, and if it is how to fix it-

I am sorry if this is a mistake, but Iw as expecting the results of the jacobians to match


#26

Hi Chirs -
Sorry to bother you once again -
I wrote a little symbolic function to calculate jacobians (see example) - Even though it calculates them correctly (I used one of my prior example problems to have direct comparison) When I pass the jacobian function to my ODE model it crashes - however it does not without the jacobian.
I am lost, be trying to figure out the issue for hours - and I could use some advice if you have time.

as always thanks and here is the code!


using Pkg
using DifferentialEquations
using SymPy
using Sundials
using LinearAlgebra
using Plots; plotly()


p=[-2,1.2,-1.0]
dIC=[-3.2,0]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]


function jacobi2(p, varnum)
    U = [symbols("u$i") for i in 1:varnum]
    A = [p[1]  -p[2]*U[1] 0;
     U[2]  p[3] 0;
    1 1 -1]
    J=jacobian(A*U,U)
    Jc= lambdify(J,U)
    #J1=N(J)
    Jc
end


function jacobi3(J,u,p,t)
    J=jacobi2(p,3)(u...)
    J=convert(Array{Float64,2}, J)
#     print(typeof(J))
    print(J)
    nothing
end

M = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]





function f5(du,u,p,t)
     
    
    A = [p[1] -p[2]*u[1] 0;
     u[2]  p[3] 0
    1 1 -1]
    mul!(du,A,u) 
    
end



test = ODEFunction(f5;mass_matrix=M,jac=jacobi3)
testprob = ODEProblem(test,IC,tspan,p)
soltest = DifferentialEquations.solve(testprob,Rodas5())

SingularException(3)

Stacktrace:
[1] checknonsingular at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\factorization.jl:12 [inlined]
[2] #lu!#103(::Bool, ::Function, ::Array{Float64,2}, ::Val{true}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lu.jl:41
[3] lu! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lu.jl:37 [inlined] (repeats 2 times)
[4] (::LinSolveFactorize{typeof(lu!)})(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Bool) at C:\Users\awolf-yadlin.julia\packages\DiffEqBase\ZQVwI\src\linear_nonlinear.jl:9
[5] perform_step!(::OrdinaryDiffEq.ODEIntegrator{Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType},true,Array{Float64,1},Float64,Array{Float64,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType},OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},LinSolveFactorize{typeof(lu!)},Nothing,ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Float64,1},1}}}}},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},LinSolveFactorize{typeof(lu!)},Nothing,ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Float64,1},1}}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}, ::OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},LinSolveFactorize{typeof(lu!)},Nothing,ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Float64,1},1}}}, ::Bool) at C:\Users\awolf-yadlin.julia\packages\OrdinaryDiffEq\miOSH\src\perform_step\rosenbrock_perform_step.jl:1176
[6] perform_step! at C:\Users\awolf-yadlin.julia\packages\OrdinaryDiffEq\miOSH\src\perform_step\rosenbrock_perform_step.jl:1124 [inlined]
[7] solve!(::OrdinaryDiffEq.ODEIntegrator{Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType},true,Array{Float64,1},Float64,Array{Float64,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType},OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},LinSolveFactorize{typeof(lu!)},Nothing,ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Float64,1},1}}}}},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},LinSolveFactorize{typeof(lu!)},Nothing,ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},Float64},Float64,1},1}}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}) at C:\Users\awolf-yadlin.julia\packages\OrdinaryDiffEq\miOSH\src\solve.jl:345
[8] #__solve#257(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\awolf-yadlin.julia\packages\OrdinaryDiffEq\miOSH\src\solve.jl:7
[9] __solve at C:\Users\awolf-yadlin.julia\packages\OrdinaryDiffEq\miOSH\src\solve.jl:6 [inlined] (repeats 3 times)
[10] #solve#429(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType}) at C:\Users\awolf-yadlin.julia\packages\DiffEqBase\ZQVwI\src\solve.jl:39
[11] solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f5),Array{Float64,2},Nothing,Nothing,typeof(jacobi3),Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Rodas5{0,true,LinSolveFactorize{typeof(lu!)},DataType}) at C:\Users\awolf-yadlin.julia\packages\DiffEqBase\ZQVwI\src\solve.jl:27
[12] top-level scope at In[31]:56

You can see from the snippet below that jacobi2 returns the right values given a vector u - so I am confused as why would this happen? I assume the issue is somewhere in the conversion from symbolic to numeric, but in all honesty I am at a loss- The behavior also occurs when using DAEs.

ex= convert(Array{Float64,2},jacobi2(p,3)(IC...)) returns 

3×3 Array{Float64,2}:
-3.2 -1.2 0.0
1.0 0.0 0.0
1.0 1.0 -1.0

Now with DAEs


function jacobi4(J,du,u,p,gamma,t)
    J=jacobi2(p,3)(u...)
    J=convert(Array{Float64,2}, J)
    J[1,1] = gamma - J[1,1]
    J[2,2] = gamma - J[2,2]
    println(J)
    nothing
end

fg2 = DAEFunction(f1;jac=jacobi4)

dIC=[-3.2,0]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]

testProb2 = DAEProblem(fg2,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = Sundials.solve(testProb2,IDA())

Return:

[3694.55 -1.19896 0.0; 1.0 3691.35 0.0; 1.0 1.0 -1.0] [14768.6 -1.19974 0.0; 1.0 14765.4 0.0; 1.0 1.0 -1.0] [59064.8 -1.19993 0.0; 1.0 59061.6 0.0; 1.0 1.0 -1.0] [2.3625e5 -1.19998 0.0; 1.0 2.36246e5 0.0; 1.0 1.0 -1.0] [9.44989e5 -1.2 0.0; 1.0 9.44986e5 0.0; 1.0 1.0 -1.0] [3.77995e6 -1.2 0.0; 1.0 3.77994e6 0.0; 1.0 1.0 -1.0] [1.51198e7 -1.2 0.0; 1.0 1.51198e7 0.0; 1.0 1.0 -1.0] [6.04791e7 -1.2 0.0; 1.0 6.04791e7 0.0; 1.0 1.0 -1.0] [2.41916e8 -1.2 0.0; 1.0 2.41916e8 0.0; 1.0 1.0 -1.0] [9.67665e8 -1.2 0.0; 1.0 9.67665e8 0.0; 1.0 1.0 -1.0]

[IDAS ERROR] IDASolve At t = 0 and h = 2.58354e-010, the corrector convergence failed repeatedly or with |h| = hmin.

Out[27]:

retcode: ConvergenceFailure Interpolation: 3rd order Hermite t: 1-element Array{Float64,1}: 0.0 u: 1-element Array{Array{Float64,1},1}: [1.0, 1.0, 2.0]


#27

Forgot the DAE equations and Real jacobian to complete post 26

function f1(res,du,u,p,t)
res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
res[3]= u[2] + u[1] - u[3]
end

function g1(J,du,u,p,gamma,t)
J[1,1] = gamma - p[1] + p[2] * u[2]
J[1,2] = p[2]* u[1]
J[1,3] = 0
J[2,1] = - 1 * u[2]
J[2,2] = gamma - p[3] - u[1]
J[2,3] = 0
J[3,1] = 1
J[3,2] = 1
J[3,3] = -1
print(J)
nothing
end

fg1 = DAEFunction(f1;jac=g1)

dIC=[0.8,2]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]

testProb2 = DAEProblem(fg1,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = Sundials.solve(testProb2,IDA())


#28

I ran

using Sundials

function f1(res,du,u,p,t)
    res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
    res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
    res[3]= u[2] + u[1] - u[3]
end

function g1(J,du,u,p,gamma,t)
    J[1,1] = gamma - p[1] + p[2] * u[2]
    J[1,2] = p[2]* u[1]
    J[1,3] = 0
    J[2,1] = - 1 * u[2]
    J[2,2] = gamma - p[3] - u[1]
    J[2,3] = 0
    J[3,1] = 1
    J[3,2] = 1
    J[3,3] = -1
    print(J)
    nothing
end

fg1 = DAEFunction(f1;jac=g1)

p=[-2,1.2,-1.0]
dIC=[0.8,2]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]

testProb2 = DAEProblem(fg1,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = Sundials.solve(testProb2,IDA())

and it seemed to work fine?


#29

Yes sorry that is an example of the function working-
The issue is when I solve for the jacobian symbollically (post 26 I think)

using SymPy
using Sundials

function f1(res,du,u,p,t)
res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
res[3]= u[2] + u[1] - u[3]
end

function g1(J,du,u,p,gamma,t)
J[1,1] = gamma - p[1] + p[2] * u[2]
J[1,2] = p[2]* u[1]
J[1,3] = 0
J[2,1] = - 1 * u[2]
J[2,2] = gamma - p[3] - u[1]
J[2,3] = 0
J[3,1] = 1
J[3,2] = 1
J[3,3] = -1
print(J)
nothing
end

fg1 = DAEFunction(f1;jac=g1)

dIC=[0.8,2]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]


# The idea here is to have a function that symbollicaly calculates the jacobian and then # pass that jacobian to the DAE or ODE problem - but something in the type I assume # is not working.

function jacobi2(p, varnum)
    U = [symbols("u$i") for i in 1:varnum]
    A = [p[1]  -p[2]*U[1] 0;
     U[2]  p[3] 0;
    1 1 -1]
    J=jacobian(A*U,U)
    Jc= lambdify(J,U)
    #J1=N(J)
    Jc
end


function jacobi4(J,du,u,p,gamma,t)
    J=jacobi2(p,3)(u...)
    J=convert(Array{Float64,2}, J)
    J[1,1] = gamma - J[1,1]
    J[2,2] = gamma - J[2,2]
    nothing
end

fg2 = DAEFunction(f1;jac=jacobi4)

dIC=[-3.2,0]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]

# Works -> Explicit Jacobian
testProb1 = DAEProblem(fg1,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = Sundials.solve(testProb1,IDA()) #works well 

# Symbolically Calculated Jacobian don't work
testProb2 = DAEProblem(fg2,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = Sundials.solve(testProb2,IDA()) # does not work well

# Here is just to show that jacobi2 works -> it is properly evaluated at IC

convert(Array{Float64,2},jacobi2(p,3)(IC...))

So I am not sure why jacobi4 or jacobi3 in post 26 (when using ODE instead of DAE fails)

Hope this clarifies my confusion (sounds like an oxymoron)

A


#30

I do think I made some progress -
But I am still having convergence problems for the DAE and I can not figute why as my symbolic jacobias is the same as in function g1 above -


using SymPy
using Sundials
using DifferentialEquations

function f1(res,du,u,p,t)
res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
res[3]= u[2] + u[1] - u[3]
end

# Equations
#du[1]  =  p[1] * u[1] - p[2] * u[1]*u[2]
#du[2]  = p[3] * u[2] + u[1]*u[2]
#0  =  -u[2] - u[1] + u[3]




# and jacobian is calulcated as: gamma*df1/d(du) + df1/du
# whihc is the same as  gamma*df1/d(du) - J(Equations), whihc is what I
# calculate using sympy

function g1(J,du,u,p,gamma,t)

# only included here for completion - I am using jacobi5 as jacobian function
J[1,1] = gamma - p[1] + p[2] * u[2]
J[1,2] = p[2]* u[1]
J[1,3] = 0
J[2,1] = - 1 * u[2]
J[2,2] = gamma - p[3] - u[1]
J[2,3] = 0
J[3,1] = 1
J[3,2] = 1
J[3,3] = -1
print(J)
nothing
end


p=[-2,1.2,-1.0]
dIC=[0.8,2]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]


# The idea here is to have a function that symbollicaly calculates the jacobian and then # pass that jacobian to the DAE or ODE problem - but something in the type I assume # is not working.

varnum=3
U = [symbols("u[$i]") for i in 1:varnum]
A = [p[1]  -p[2]*U[1] 0;
 U[2]  p[3] 0;
1 1 -1]
J1=jacobian(A*U,U)
Jc=lambdify(J1,U)

    

function jacobi5(J,du,u,p,gamma,t)
    J=Jc(u...)
    J=convert(Array{Float64,2}, J)
    # as per def of jacoboian for DAE
    #
    J= -J
    J[1,1] = gamma  + J[1,1]
    J[2,2] = gamma + J[2,2]
    nothing
end

fg3 = DAEFunction(f1;jac=jacobi5)
testProb3 = DAEProblem(fg3,dIC,IC,tspan,p,differential_vars=diff_vars)
sol3= Sundials.solve(testProb3, IDA())

When I run this no more awful errors arise,

[IDAS ERROR] IDASolve
At t = 0 and h = 2.58354e-10, the corrector convergence failed repeatedly or with |h| = hmin.

or

[IDAS ERROR] IDASolve
The value tstop = 10 is behind current t = 0, in the direction of integration.

and I am not sure why as the jacobian form Jacobi5 is the same expression as g1 which convrges


#31

This function isn’t changing the Jacobian at all. You’re creating a new array instead of modifying the existing one.


#32

So would you just then looo over j and assign J[i,j] =Jc(u…)[i,j]

I assume that would update the Jacobian instead of creating new one correct?

Thanks!


#33

Got it to work!!!
That clue about updating the jacobian made all the difference!

Here is the code in case it’s useful to you!

using SymPy
using Sundials
using DifferentialEquations

function f1(res,du,u,p,t)
res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
res[3]= u[2] + u[1] - u[3]
end

# Equations
#du[1]  =  p[1] * u[1] - p[2] * u[1]*u[2]
#du[2]  = p[3] * u[2] + u[1]*u[2]
#0  =  -u[2] - u[1] + u[3]




# and jacobian is calulcated as: gamma*df1/d(du) + df1/du
# whihc is the same as  gamma*df1/d(du) - J(Equations), whihc is what I
# calculate using sympy

function g1(J,du,u,p,gamma,t)

# only included here for completion - I am using jacobi5 as jacobian function
J[1,1] = gamma - p[1] + p[2] * u[2]
J[1,2] = p[2]* u[1]
J[1,3] = 0
J[2,1] = - 1 * u[2]
J[2,2] = gamma - p[3] - u[1]
J[2,3] = 0
J[3,1] = 1
J[3,2] = 1
J[3,3] = -1
nothing
end


p=[-2,1.2,-1.0]
dIC=[0.8,2]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]


# The idea here is to have a function that symbollicaly calculates the jacobian and then # pass that jacobian to the DAE or ODE problem - but something in the type I assume # is not working.

varnum=3
U = [symbols("u[$i]") for i in 1:varnum]
A = [p[1]  -p[2]*U[1] 0;
 U[2]  p[3] 0;
-1 -1 1]
J1=jacobian(A*U,U)
Jc=lambdify(J1,U)

dv=sum(diff_vars)
s=size(Jc(IC...))
function jacobi5(J,du,u,p,gamma,t)
    for i in 1:s[1]
        for j in 1:s[2]
            if (j== i & i≤dv)
                J[i,j]=gamma + 1*float(-1*Jc(u...)[i,j])
            else
                J[i,j]=1*float(-1*Jc(u...)[i,j])
            end
        end
    end
end

fg3 = DAEFunction(f1;jac=jacobi5)
testProb3 = DAEProblem(fg3,dIC,IC,tspan,p,differential_vars=diff_vars)
sol3= Sundials.solve(testProb3, IDA())

fg1 = DAEFunction(f1;jac=g1)
testProb1 = DAEProblem(fg1,dIC,IC,tspan,p,differential_vars=diff_vars)
sol1= Sundials.solve(testProb1, IDA())


using Plots
plot(sol1, layout=(3,1), label=["g1 X1" "g1 X2" "g1 X3"], lw=3)
plot!(sol3,linestyle = :dash, layout=(3,1), label=["Jacobi5 X1" "Jacobi5 X2" "Jacobi5 X3"], lw=3)

thanks a ton! -
A


#34

Hi -
I hope you are doing well -
If you have time, I have one last question - I am trying to ensure that all values in my solution vector are >=0. I have tried teh positive domain callback as well as as “isoutofdomain=(u,p,t) -> any(x -> x < 0, u)” solution without mayor success -
Do you have any advice on how to use these tools?
I need imtegration to stop when withe u[16] or u[19] reach zero as I ahve some unctions that depends in their square root.

Here is the code I have tried -
Thanks fro your help!

condition(u,t,integrator) = t in ustrip.([14u"d",28u"d",42u"d"])


function affect!(integrator)
    integrator.u[1] += ustrip(Dose)
    integrator.u[11] = 4#(4*ustrip(Dose))+u[11]*(u[1]+u[2]+(u[3]+u[4]+u[5])*u[16]))/(ustrip(Dose)+(u[1]+u[2]+(u[3]+u[4]+u[5])*u[16]))
end

condition1(u,t,integrator) = u[16] < 1e-8
#condition1(u,t,integrator) = u[16] 
affect1!(integrator) = terminate!(integrator)

cb = DiscreteCallback(condition,affect!)
#cb1 = ContinuousCallback(condition1, affect1!, rootfind = true)
cb1 = DiscreteCallback(condition1, affect1!)
#cb2=PositiveDomain()
cbset = CallbackSet(cb1,cb)
#cbset=CallbackSet(cb2,cb)


ODE_test_cb_sol = 
DifferentialEquations.solve(ODE_test_Prob,Rodas5(),callback=cbset,tstops=ustrip.([14u"d",28u"d",42u"d"]))

# DifferentialEquations.solve(ODE_test_Prob,Rodas5(),isoutofdomain=(u,p,t) -> any(x -> x < 0, u),callback=cb,tstops=ustrip.([14u"d",28u"d",42u"d"]))

As you can see I have tried discrete and continuous callbacks, as well as PositiveDomain and isoutofdomain -
I am not sure if I am using these tools correctly - I thought I was
Please see below the error I usually get -The code above terminating integration when U[16] <1e-8 is the only solution that had worked for me, but ideally I’d like to get to a solution where integartion stops exactly at 0 or closer than1e-8

Thanks

DomainError with -1.7502393112149852e-6:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.


#35

You need to make your derivative function compatible with negative values, and then isoutofdomain can work. Right now it’s failing in your function before it’s able to check if it’s out of the domain


#36

Ok so you mean my problem is in the definition of my equations? for example eqs 19-21?
How do you make this compatible with negative values do you use abs or multiply the equation by -1?
sorry I am a little confused by as how to proceed given that u21 will always depend on the square root of u19-

Thanks

res[19]=-62035.04908994*u[16]^(1/3) + u[19]

res[20] =-6*P[8]/u[19]^2 + u[20]

res[21]= -6*P[21]/u[19]^2 + u[21]


#37

Yes, use abs or something so it doesn’t error, and instead continues and rejects the step.


#38

Ok
thanks-

Ale


#39

Hi Chris, (I assume is likely you answering),
I have completed my model and it is know stable - I have 21 equations (15 are differential equations and 6 algebraic) - I am using ODEfunction with mass matrix, passing jac and if needed paramjac -
I am interested in starting sensitivity analysis, but I had a couple of problems with it - Following your tutorial in the ODE page this is what I did -


# this solves no problem
ODE_test = ODEFunction(ODE_Sym;mass_matrix=MM,jac=Jac_ODE_Sym, paramjac=paramJac_ODE_Sym)
ODE_test_Prob = ODEProblem(ODE_test,IC,tspan,p)
ODE_test_sol = DifferentialEquations.solve(ODE_test_Prob,Rodas5())

# this crashes with a huge error I forgot to copy sorry, but it starts like this
# MethodError: no method matching ##2379(::Float64, ::Float64, ::Float64, ::Float64, ::Float64, 
# ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, 
# etc

ODE_test_Prob_sensitivity = ODELocalSensitivityProblem(ODE_test,IC,tspan,p)
ODE_test_sol_sensitivity = DifferentialEquations.solve(ODE_test_Prob_sensitivity,Rodas5())


# this declared the system unstable after one iteration *length(u) =1  
# ODE_test_sol_sensitivity[1] = 693-element Array{Float64,1}
#  Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true # # # 
# solution is unstable. @ DiffEqBase # 
# C:..\.julia\packages\DiffEqBase\ZQVwI\src\integrator_interface.jl:156

ODE_test_Prob_sensitivity = ODELocalSensitivityProblem(ODE_test,IC,tspan,p)
ODE_test_sol_sensitivity = DifferentialEquations.solve(ODE_test_Prob_sensitivity,DP8())


# this been running for 1.5 hrs now
ODE_test_Prob_sensitivity = ODELocalSensitivityProblem(ODE_test,IC,tspan,p)
ODE_test_sol_sensitivity = DifferentialEquations.solve(ODE_test_Prob_sensitivity,Vern9())

# update vern9() returned this error
# DomainError with -0.00018137762147211046:
# log will only return a complex result if called with a complex argument. Try log(Complex(x)). 
# so i guess need to ensure all params are real?

I am just curious is this expected behavior? Should the same solver Rodas5() that solved the system no problem whenusing ODEProblem be expected to solve also when running ODELocalSensitivityProbelm analysis?

Maybe I am doing something wrong when defining the ODELocalSensitivityProbelm? Or is the issue a function of having algebraic equations in the system? If so is there a way you recomend to address that issue?

Would this or the Morris/sobol methods benefit from using jac_param explicitly when defining ODEFunction?

thanks a ton!

Ale


#40

I don’t think this is what you want to be using. Why not use ForwardDiff or adjoint methods? For non-trivial mass matrices this might not be correct anyways, we should make it throw an error in that case. Please open an issue on that so we can remember to throw an error.