Hello, I am a R user, and just start to learn Julia. I used to use R ‘desolve’ package to solve systems of ODEs, but the speed of Julia ‘DifferentialEquations’ make me switch to Julia to handle some large systems.
My questions is how to estimate the jacobian matrix of the following ODE systems at end solutions, which is usefull to infer local stability of this system. It is quite straightfoward to do this in R by ‘rootSolve::jacobian.full’. I am wondering if there is a similar function to do it easily in Julia. I found funcitons like ‘ForwardDiff.jacobian’, and some discussion on ODEs like LORENZ, but I did not know how to use it on my ODEproblem. I am quite unfamiliar with Julia script. I will greatly appreciate it if some one can provide a testable example for my following ode problem. Thank you very much.
#my code:
using Random, Distributions, LinearAlgebra,DifferentialEquations,Plots
p_p_m = rand(Normal(0,0.1), 10,10)
p_p_m[diagind(p_p_m)].=-rand(Uniform(0.5,1.5),10)
m_m_m = rand(Normal(0,0.1), 10,10)
m_m_m[diagind(m_m_m)].=-rand(Uniform(0.5,1.5),10)
p_m_m=abs.(rand(Normal(0,0.5), 10,10))
m_p_m=abs.(rand(Normal(0,0.5), 10,10))
p_R=rand(Uniform(0.5,1.5),10)
m_R=rand(Uniform(0.5,1.5),10)
p_i=1:10
m_i=11:20
parms=Array[p_R,m_R,p_p_m,m_m_m,p_m_m,m_p_m,p_i,m_i]
function pm_model!(du,u,parms,t)
p=u[p_i]
m=u[m_i]
du[p_i] = p.*( parms[1] +parms[3]*p+ parms[5]m./(0.5parms[5]m.+1))
du[m_i] = m. ( parms[2] +parms[4]*m+ parms[6]p./(0.5parms[6]*p.+1))
end
u0 = rand(Uniform(0.5,1.5),20)
tspan = (0.0,1000.0)
prob = ODEProblem(pm_model!,u0,tspan,parms)
sol = solve(prob)
plot(sol)