thanks
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
Look at http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1
Thanks again
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(p) 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
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]
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())
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?
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
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
This function isnāt changing the Jacobian at all. Youāre creating a new array instead of modifying the existing one.
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!
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
1 Like
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.
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
1 Like
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]
Yes, use abs or something so it doesnāt error, and instead continues and rejects the step.
Ok
thanks-
Ale
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
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.