Hi. I have a system of differential equations containing 62500 simultaneous equations. The system dynamics are linear and of the form: [dU/dt] = [Jacobian]x[U].
The jacobian is sparse (non-zero vals ~ 0.6%) and banded. However, the system is also highly stiff. The largest value in the Jacobian matrix is roughly 1E12 and the smallest nonzero value is probably around 1E3 so I’m using an implicit solver: CVODE_BDF with GMRES as the linear solver. However, this method results in convergence failures from the solver. I have tried reducing the tolerances till reltol =1E-10 and abstol = 1E-7. I’m still not obtaining satisfactory results and convergence failures still occur. Currently, I am trying alternative solvers such as: QNDF, FBDF etc., . I will update on those results soon.
Note: This problem is a reduced version of a larger original problem. The original problem contains ~ 500k - 1 Million simultaneous differential equations of similar form. Also, I dont require the full solution at all points in time. I am only interested in the total amount of species with time: i.e. sum(u[:]) with respect to t.
A brief description of the physics: I am trying to solve a nanometric species transport problem. Assume a nanometric cuboid box (25nm x 25nm x100nm), where the bottom face of the box has an initial distribution of species. With time, the species hop and transport to the top face of the box. Once the species reach the top face, they exit (the exit is modelled as decay when they reach the top face).
Do you need the full solution? If you only need the endpoint this is just an expv from ExponentialUtilities.jl call to compute the solution at an arbitrary time. Even if that’s not the case, you might want to try an exponential integrator such as ODE Solvers · DifferentialEquations.jl.
No, I dont require the full solution (i.e. solution of indvidual differential equations at all points in time). I am only interested in the total amount of species with time
i.e. sum(u[:]) with respect to t.
Also, are exponential integrators a good approach? I checked the documentation, it says they only operate on a fixed step-size (Please correct me if im wrong). The fastest dynamics in my system is 1E12 Hz and the slowest is about 1E3 Hz. So would’nt the solver struggle to march forward in time? Kindly clarify.
They work on a fixed time step, but it is not an explicit solver where your time step should be around 1e-12 to be stable. Exponential integrator are based on the idea of solving exactly for the linear part of the ODE (in this case this applies to the full system since you said it is linear). This means writing the solution of your ODE
du/dt = Au
as
u(t) = e^{tA}u_0
where u_0 is some initial condition. So really, in your case, it is enough that you compute the matrix exponential for the time steps of interest to compute the solution. Exponential integrators basically do this but when there are other nonlinear terms they are treated in a similar way to other numerical methods but after solving for the linear part with the exponential of a matrix.
I knew that this analytical solution exp(Jacobian*t)U0 exists. I was’nt sure if i had to solve the problem by taking this into the Eigen-coordinate space and representing the decoupled dynamics as z(t) = exp(lamdat)*z0. However, for a large jacobian matrix such as in this case the computation of the Eigen-values and vectors is highly impractical. Thats where my confusion was.
Thanks for clarifying. I will try out exponential integrators and let you know the results.
I have one more doubt as well The matrix exponential is just the taylor’s series expansion of the jacobian matrix, so you will have quadratic, cubic, quartic etc… terms which may not be sparse even though the jacobian is sparse. So wouldnt the RAM be chocked because the solver has to store a dense exp(At) matrix in RAM to solve the problem. Kindly clarify this too.
The exponential integrator will basically project the jacobian on a well chosen basis (function of u0) and return u(t) at specified times. It is designed to work in very large dimensions.
The matrix exponential is rarely computed using the Taylor Series due to the reasons that you are outlining. If you check the previous library suggested to you in one of the replies to your post, you can see for yourself the implementation on how to compute these matrix exponentials.
The piece that you are missing is that you don’t actually need to compute exp(J) to compute exp(J)*u0. exp(J) as you mentioned is a dense matrix, but you can compute the action of exp(J) on a vector with only matrix vector multiplies. Note that the following is not the “good” way to do this (the good ways use things like krylov subspaces and other fancyness), but it should make it clear how the density problem is avoided.
function expv(A, u0)
out = zeros(length(u0))
tmp = copy(u0)
i = 1
while true
tmp .= A*tmp ./ i
out .+= tmp
i += 1
if norm(tmp) < 1e-8*norm(out) # if converged
return out
end
end
end
Basically, rather than computing (A + A^2/2 +A^3/6...)*u0, we compute A*u0 + A(A*u0)/2+A(A(A*u0))/6... which ensures that we never need to explicitly store the dense matrix.
RESULT: expv gives correct results only until a t_end value of 1e-9 seconds, beyond which the values are 0 everywhere, but i need results till a time-span of atleast 1e-3 seconds.
Note: The max values in the jacobian are of the order 1e12, the minimum nonzero values are 1e-3 and i want to solve the problem till a span of 1e-3 seconds (a large time-span considering the fastest dynamics are 9 orders faster than the slowest dynamics)
RESULT: expv_timestep gives more accurate solutions but the ‘dt’ value was not adaptive (although i set adaptive= true) and the ‘dt’ value was always around 1e-10 seconds (I simulated for about 2-3 days and i only got until 1e-6 seconds. CVODE_BDF did much better and was able to solve until 200e-6 seconds within 20 hours of simulation time. However CVODE_BDF too became unstable later on)
Will update the results of adaptive exponential RK methods exprb32 and exprb43 soon.
Also, I did not want to try fixed time step exponential methods such as exponential RK and exponential propagation iterative methods because I want the results over a large timespan.
Hi, thanks a lot for helping me clear this problem. As i mentioned before,: I am trying to solve a nanometric species transport problem. Where the species have to travel from the bottom face of a cuboid box (25nm x 25nm x100nm) to the top face (where they exit). I inspected the jacobian matrix carefully and i came to understand that at some elements, the influx rate was orders higher than the outflux rate. This did not let the possibility of species decay/transport out of that element (if they flow out they will be pushed back into that element). This probably resulted in undamped oscillations in the system’s dynamics. This was the main reason for ODE solver failure.
I Identified the specific locations where this occured and added an additional decay mode to the existing dynamics to force the species to decay in such elements. Now after doing the same, CVODE_BDF works just fine. QNDF also works fine. Exponential solvers Exprb32 and 43 are extremely slow for the working problem. It was more of a Physics related issue than an ODE-related issue (Sorry for being the Patrick)
Also, I work on a lot of physics-based modelling for high voltage applications, FEM and wireless power transfer. So if you had any problems in such domains, I would be happy to reciprocate.
The matrix exponential is just the taylor’s series expansion of the jacobian matrix, so you will have quadratic, cubic, quartic etc… terms which may not be sparse even though the jacobian is sparse. So wouldnt the RAM be chocked because the solver has to store a dense exp(At) matrix in RAM to solve the problem. Kindly clarify this too.