rjhogan / Adept-2

Combined array and automatic differentiation library in C++
http://www.met.reading.ac.uk/clouds/adept/
Apache License 2.0
163 stars 29 forks source link

Question about the Adept-2 #23

Closed jianqixi closed 1 year ago

jianqixi commented 2 years ago

Dear All,

I am trying to build the Jacobian matrix for one ODE question, like dy/dt=f(t,y). I have programmed the C++ code for the ODE functions. And I want to build one Jacobian matrix of df/dy for the ODE solver. I know I can create the hand code of df/dy by myself, but it is not effective. So I want to know whether the Adept algorithm can help me to build this Jacobian matrix? The size of this matrix is around 3000*3000. Thank you,

rjhogan commented 2 years ago

Yes you could use Adept for that - you would need to code up a function that computes f(t,y) using Adept's active types, e.g. adept::adouble, adept::aVector etc. Have a look at the two listings in section 2.4 of the documentation (http://www.met.reading.ac.uk/clouds/adept/adept_documentation_2.1.pdf) - one uses Adept's scalar types (adouble), the other its array types.

jianqixi commented 2 years ago

Thank you for your quick reply. As a beginner of this kind of algorithm, I am wondering do you have some specific examples because I am reading the user guide, but it looks like it has multiple examples but maybe not be related to mine. So I wonder would you like to give me some specific examples. Thank you very much!

rjhogan commented 2 years ago

The test directory of the package is full of examples. The function beginning on line 117 here: https://github.com/rjhogan/Adept-2/blob/master/test/test_minimizer.cpp#L117 computes a Jacobian matrix of a function of the form y(x), where y and x are vectors. This seems to be almost identical to your problem. Obviously my program does a bit more with the result - in fact this code was used to produce the data for the animation here: http://www.met.reading.ac.uk/clouds/adept/optimization.html

jianqixi commented 2 years ago

Thank you for your recommendations.

jianqixi commented 2 years ago

@rjhogan I have another question about the usage of Adept. From the manual, it seems like the AD algorithm in Adept will provide the numerical value for each element in the Jacobian matrix. But if there is a variable in the element, could I also use this algorithm? One specific example is that:

One of the ODE equations in my case is dy(i)/dt = Ay(i) + By(i-1)y(i+1), then the element in the Jacobian matrix could be By(i+1), but if the y(i+1) is the unknown which would be the solution of from the solver. In this case, could I still use the Adept numerical AD algorithm? Thank you for your clarification.

rjhogan commented 2 years ago

No you can't use Adept for that - you would need a symbolic math language. I wonder if Adept is the tool you need - Adept is for differentiating very complex functions, whereas you know the expression for the derivative of a function and you want to integrate it to work out the function at a new time. If you want to evolve y(i) in time then the classic approximate approach would be to chose a timestep and use something like the Midpoint Method or (for more accuracy) Runge-Kutta-4. This is useful if the right-hand-side is too complicated - in your case it might be simple enough that there is an analytic solution but I don't know how to handle the nonlinear "B" term (I should say I'm a physicist not an applied mathematician). If the right hand side consisted only of linear terms [nothing with y(i)*y(j)] then there is an analytic solution in terms of matrix exponentials.

jianqixi commented 2 years ago

Thank you for your clarification! Yes, I am currently trying to use a linear matrix solver (Dense Matrix solver) for the ODE function, f(y,t). The current solver definitely has a default Jacobian matrix by calling the ODE function multiple times. But using the default method is very time-consuming since my matrix size is too large (~2886).

I agree with you that the nonlinear term in the right-hand may not easy to be solved analytically. The current solver that I used is based on the BDF algorithm in the Newton matrix, it could be based on the numerical method to solve the nonlinear issue. But to be honest, I also do not fully understand how it works.

rjhogan commented 2 years ago

BDF is solving the same family of problems as Runge-Kutta but I didn't know that (unlike Runge-Kutta) BDF makes use of df/dy - I can see how it might help extrapolate to the new value of t. Let me know if you need help using Adept to compute df/dy. I guess not all of your ODEs are as simple as f(i)=Ay(i) + By(i-1)y(i+1) because the Jacobian of this is trivial to code up by hand: df(i)/d[y(i-1), y(i), y(i+1)] = [By(i+1), A, By(i-1)] where I'm using square brackets to indicate a vector.

jianqixi commented 2 years ago

Really appreciated your help. I wonder when using Adept for the calculation of df/dy, is it symbolic? If it is, that would be very helpful. Thank you!

rjhogan commented 2 years ago

No, it is algorithmic differentiation, as distinct from symbolic differentiation and numerical differentiation - see https://en.wikipedia.org/wiki/Automatic_differentiation