Nlsolve for find root for a system of ODE for an implicit runge kutta method

This is an ODE solver, The t, dt, u and the length of k changes based on the ODE function we are trying to solve. In this particular example (van der Pol equation) there are 2 variables: u[2, 0] at t =0. Therefore k takes 2 x 3 ones (n*3) …n: number of variables, 3 is the number of unknown for the Radau IIA 5th order solver.

But if I change the ODE function to that of for example the Robertson Chemical Reactions, we have 3 variables: u =[1, 0, 0] and therefore k becomes 3 x 3 ones: k=[1,1,1,1,1,1,1,1,1]

The f function is modifying the ODE function fu to accommodate the Radau IIA 5th Order need for 3 unknowns in each step; there are 3 steps to the Radau5 solver.

See this article on the Radau IIA 5th order butcher tableau: wiki

Therefore I have set up my code as follows:

ODE function for van der Pol equation:

function fu(u, t)
  p = 1000
  du = [
    u[2],
    p * (1 - u[1]^2) * u[2] - u[1]]
  return du
end

The initial conditions (t, u, dt):

t = 0.0 # begin
dt = 0.001 # timestep
tn = 3000 # end

u = [2.0, 0.0]

And then the f function to modify the fu function to accommodate the Radau5 3 steps with each step having 3 unknowns:

function f(u, t, dt, k)
  n = length(u)

  a00 = 11.0/45.0 -  (7.0 * sqrt(6.0)) / 360.0
  a01 = 37.0/225.0 -  (169.0 * sqrt(6.0)) / 1800.0
  a02 = -2.0/225.0 +  sqrt(6.0) / 75.0
  a10 = 37.0/225.0 +  (169.0 * sqrt(6.0)) / 1800.0
  a11 = 11.0/45.0 +  (7.0 * sqrt(6.0)) / 360.0
  a12 = -2.0/225.0 -  sqrt(6.0) / 75.0
  a20 = 4.0/9.0 -  sqrt(6.0) / 36.0
  a21 = 4.0/9.0 +  sqrt(6.0) / 36.0
  a22 = 1.0/9.0

  c0 = 2.0/5.0 - sqrt(6.0) / 10.0
  c1 = 2.0/5.0 + sqrt(6.0) / 10.0
  c2 = 1.0

  
  out = ones(n*3)
  xcurr = zeros(n)
    
  for i in 1:n
    j = (i + (i * 2)) - 2
    xcurr[i] = u[i] + (a00 * k[j] + a01 * k[j+1] + a02 * k[j+2])
  end

  o = fu(xcurr, t + c0 * dt)
  for i in 1:n
    j = (i + (i * 2)) - 2
    out[j] = dt * o[i] - k[j]  
    xcurr[i] = u[i] + (a10 * k[j] + a11 * k[j+1] + a12 * k[j+2])
  end
    
  o = fu(xcurr, t + c1 * dt)
  for i in 1:n
    j = (i + (i * 2)) -2
    out[j+1] = dt * o[i] - k[j+1]  
    xcurr[i] = u[i] + (a20 * k[j] + a21 * k[j+1] + a22 * k[j+2])
  end
    
  o = fu(xcurr, t + c2 * dt)
  for i in 1:n
    j = (i + (i * 2)) -2
    out[j+2] = dt * o[i] - k[j+2]
  end

  return out
end

Now I want to solve it without having to code my own Newton’s method (which I already did and is currently in my code and working successfully):

I have k as a guess for the unknowns (because Radau5 has 3 steps with each step having 3 unknowns, we are using the 2 variables for each step):

m = length(u)
k = ones(3*m)

And now the solver:

#c1 = nrm(u, t, dt, k) # works perfectly fine (my own Newton's method code)
c1 = nlsolve(f, k) # which does not work

Can you guys help me make this nlsolve work please?

I used your example: but the error says:

ERROR: BoundsError: attempt to access 6-element Vector{Float64} at index [7]

And the location is in the f function at the “xcurr” on the first loop.

https://docs.julialang.org/en/v1/base/punctuation/
https://juliateachingctu.github.io/Julia-for-Optimization-and-Learning/dev/lecture_02/arrays/

1 Like

How does it fair vs DifferentialEquations.jl?

See, it’s still missing all of the standard tricks of the trade. For a FIRK, you can use the symmetry properties to evaluate the nonlinear equation in the complex plane:

That then greatly reduces the size of the Jacobian and improves SIMD. Then you can add Jacobian reuse strategies by exploiting the adaptive time stepping and rejecting steps which are failing.

So actual implicit RK code would never use the algorithm you’re using here.

The paper suggests matrices to solve the nonlinear equations. But I haven’t found an example of this implementation on the web. I constructed my Radau5 based on the usual implementation of the explicit methods which is readily available anywhere. I just figured out on my own how to include the ODE function within the Radau5 equation to solve via Newton. I’ve used this a lot.

I do think that the libraries in Python, Julia, and other languages have better algorithms for ODE solving and as such are faster, but they make use of, like you said “tricks of the trade”. They will get the same result I’m getting, albeit faster. My Radau5 usually solves everything for me under 1 minute or less. The only ODE that takes way too long is the van der Pol Equation. It’s taking 5 mins in Julia.

You seem to know a lot about one of the tricks they make use of in their algorithm. Do you have a personal implementation of the Radau5 with matrices?

The usually implementation, Hairer’s radau and radau5 implementation in Fortran, both use this trick.

Yes, in DifferentialEquations.jl. The code I linked is the Julia code that uses these tricks:

https://diffeq.sciml.ai/stable/solvers/ode_solve/#Fully-Implicit-Runge-Kutta-Methods-(FIRK)

On Van Der Pol it’s not a great algorithm

https://benchmarks.sciml.ai/html/StiffODE/VanDerPol.html

It solves in about 10^(-2) seconds, but high order Rosenbrock methods are about 4x faster on this problem even with using a tuned Jacobian reuse scheme (IMO the one that Hairer suggests is a bit wasteful on modern CPUs given the timing difference between LU and backsolves). I’d be curious if you had any more improvements to the algorithm beyond that.

1 Like

This code does not seem to be a standalone code, but may have dependencies and is very big. It looks like it is using the Radau 3rd Order as a way to measure accuracy.

Don’t you have something simple, like only for Radau5 but with the tricks (matrices)?

Also, as you mentioned Rosenbrock methods and Rodas3, there isn’t a single article on wikipedia for it. So I just don’t know how they are constructed.

I just looked at your profile. I’m very impressed with your PhD in maths and you’re a major Julia developer! I like maths and physics (cosmology). I would really appreciate it if you could provide me with a standalone Radau5 for low tolerance and a Rodas3 for higher tol please!

But something I could both code in Python and in c++ …I’m very fluent in Python, Sagemath & VBA, dabbled in Javascript and R…I used Mathematica & Maxima as a verifier for mathematics and as Julia is similar in syntax with Python…I’ve started out with Julia and c++…But I’m looking to work more in Julia

My Radau code as it stands with regard to the van der Pol function takes 1 hour in Python, 5 mins in Julia and longer than the other two in c++

By the way, I don’t see Dopri 8(7) with 13 stages listed in differential Equations.jl Why is that? article

There is another algorithm not listed also, but which I really needed. The Gauss-Jackson 8th Order for orbit propagation - article
and article2

It’s pretty hard to guess what you did when you say “I used your example.” Did you try to run the code I provided by cutting and pasting it exactly as is? Did you try to run nlsolve using my f! and your f? Or something else? We need a complete minimum working example and explicit error messages to be able to help.

Please do read Please read: make it easier to help you and follow its recommendations.

1 Like

Ok
Here is my code for nlsolve:

using NLsolve

function f(u, t, dt, k)
  n = length(u)

  a00 = 11.0/45.0 -  (7.0 * sqrt(6.0)) / 360.0
  a01 = 37.0/225.0 -  (169.0 * sqrt(6.0)) / 1800.0
  a02 = -2.0/225.0 +  sqrt(6.0) / 75.0
  a10 = 37.0/225.0 +  (169.0 * sqrt(6.0)) / 1800.0
  a11 = 11.0/45.0 +  (7.0 * sqrt(6.0)) / 360.0
  a12 = -2.0/225.0 -  sqrt(6.0) / 75.0
  a20 = 4.0/9.0 -  sqrt(6.0) / 36.0
  a21 = 4.0/9.0 +  sqrt(6.0) / 36.0
  a22 = 1.0/9.0

  c0 = 2.0/5.0 - sqrt(6.0) / 10.0
  c1 = 2.0/5.0 + sqrt(6.0) / 10.0
  c2 = 1.0

  
  out = ones(n*3)
  xcurr = zeros(n)
    
  for i in 1:n
    j = (i + (i * 2)) - 2
    xcurr[i] = u[i] + (a00 * k[j] + a01 * k[j+1] + a02 * k[j+2])
  end

  o = fu(xcurr, t + c0 * dt)
  for i in 1:n
    j = (i + (i * 2)) - 2
    out[j] = dt * o[i] - k[j]  
    xcurr[i] = u[i] + (a10 * k[j] + a11 * k[j+1] + a12 * k[j+2])
  end
    
  o = fu(xcurr, t + c1 * dt)
  for i in 1:n
    j = (i + (i * 2)) -2
    out[j+1] = dt * o[i] - k[j+1]  
    xcurr[i] = u[i] + (a20 * k[j] + a21 * k[j+1] + a22 * k[j+2])
  end
    
  o = fu(xcurr, t + c2 * dt)
  for i in 1:n
    j = (i + (i * 2)) -2
    out[j+2] = dt * o[i] - k[j+2]
  end

  return out
end


function fu(u, t)
  p = 1000
  du = [
    u[2],
    p * (1 - u[1]^2) * u[2] - u[1]]
  return du
end

t = 0.0 # begin
dt = 0.001 # timestep
tn = 3000 # end

u = [2.0, 0.0]

n = length(u)

k = ones(3*n)
  
    
function f!(fx, u)
   fx .= f(u, t, dt, k)
end

  c1 = nlsolve(f!, k)

And error is:

ERROR: BoundsError: attempt to access 6-element Vector{Float64} at index [7]

It is still not clear to me what mathematical problem you are solving. Van der Pol, ODE integrators, and nonlinear solving I understand. But I can’t understand what your f function or f! function represents, or what the search space of the root-finding is supposed to be.

nlsolve will find a root of the equation f(x) = 0 for a function f : \mathbb{R}^n \rightarrow \mathbb{R}^n. The code you’ve written gives nlsolve the function f! which maps the two-d u into the six-d output of f. Thus you’re trying to find a root of an \mathbb{R}^2 \rightarrow \mathbb{R}^6 function. That doesn’t make sense --there are more equations than unknowns.

Similarly, the second argument of nlsolve should be an initial guess for the zero of f, that is, the variable x for which f(x) = 0. In your code, f! takes a two-d input u, but you provide nlsolve with a six-d initial guess, for k.

What is the mathematical problem you’re trying to solve? What variable is nlsolve supposed to be searching over? What is the dimension of that variable? What does the k variable represent? Is the search space u or k?

Dopri8(7) is not an efficient tableau. This has been known for years. This is why it was replaced really quickly by the authors themselves with the dop853 method (which then had optimized software by Hairer) and no one really used Dopri8(7) again, except for people outside of the field for some reason… mostly confusing it with dop853 which is actually good.

We implemented 10,000 lines of tableaus, you can find that here:

and documented here:

https://docs.sciml.ai/dev/modules/DiffEqDevDocs/internals/tableaus/

constructDormandPrice8() is the method you refer to. But in no benchmark does that method do well.

In the end, the Verner efficient methods just dominate the field here though. You can see that in the paper we put out today: [2207.08135] Parallelizing Explicit and Implicit Extrapolation Methods for Ordinary Differential Equations

For explicit tableaus, the Verner methods are the really damn good ones, but they only were created in like 2017 and didn’t have software until very recently, so most older benchmarks don’t include them.

Ehh, it’s in the todo list:

https://github.com/SciML/OrdinaryDiffEq.jl/issues/331

Multi-step methods are generally underwhelming in the explicit case. I have almost never seen them win a benchmark for non-stiff ODEs no matter how much work we put into them. For the symplectic case, why not use a symplectic RK like Dynamical, Hamiltonian, and 2nd Order ODE Solvers · DifferentialEquations.jl ? Or if you really want a PCE symplectic method, the 16th order Gauss-Legendre is rather good and should outperform GJ quite handedly:

https://github.com/SciML/IRKGaussLegendre.jl

Yes, we need to do more benchmarking here. But I don’t see any good evidence for GJ which is why it’s been on the backburner.

3 Likes

The function f represents the following Runge Kutta implicit function:

k1 = h ∗ f( t + c1 ∗ h, y + (a11 ∗ k1 + a12 ∗ k2 + a13 ∗ k3))

k2 = h ∗ f( t + c2 ∗ h, y + (a21 ∗ k1 + a22 ∗ k2 + a23 ∗ k3))

k3 = h ∗ f( t + c3 ∗ h, y + (a31 ∗ k1 + a32 ∗ k2 + a33 ∗ k3))

y(n+1) = y + (b1 ∗ k1 + b2 ∗ k2 + b3 ∗ k3)

The As (Runge Kutta matrix), Bs (weights) and Cs (nodes) are the coefficients of the Radau IIA 5th Order.

The above equations for Implicit Runge Kutta method is an example for a single function of type:

f(u, t); u(t0) = y0

For coupled system of ODE, the 3 steps shown above gets multiplied by the number of variables in your system.

For a 3 variables system of ODE, you would need to solve a Jacobian 9 x 9 matrix. The function to which that Jacobian matrix is derived from is the f function in my code.

So in effect for Robertson Chemical Reactions ODE which has 3 variables:
The k1 ( for first variable)
k2 (first variable)
k3 ( first variable)

k1(second variable)
k2(second variable)
k3 (second variable)

k1(third variable)
k2(third variable)
k3(third variable)

The coefficients you have here is different from that of Hairer: link

What does “//” mean in Julia?

I do have Efficient Verner 9th Order. What’s the difference between the “Robust” and “Efficient” one?

I actually have the following:

Choose numerical method (1 - 19):
1 - Bogacki-Shampine 3rd Order
2 - Runge Kutta 4th Order
3 - Gill’s Method 4th Order
4 - Dormand Prince 5th Order
5 - Dormand Prince 8th Order
6 - Radau IIA 5th Order
7 - Adams-Bashforth-Moulton 5th Order
8 - Adams-Bashforth-Moulton 8th Order
9 - Efficient Verner 9th Order
10 - Curtis 8th Order
11 - Lobatto IIIC 4th Order
12 - Gauss-Legendre 6th Order
13 - Backward Differentiation Formula 1st Order
14 - Backward Differentiation Formula 2nd Order
15 - Backward Differentiation Formula 3rd Order
16 - Backward Differentiation Formula 4th Order
17 - Backward Differentiation Formula 5th Order
18 - Lobatto IIIA 4th Order
19 - Trapezium Rule
20 - Leapfrog Yoshida 4th Order (choose nbodyv7 or nbodyv7a)

I am using Yoshida 4th Order as referenced here and Dopri 8(7) 13 for nbody problem.

I am interested in your Yoshida 6th order.

Yes, I mentioned that Dopri8(7) is not dop853. They are separate algorithms. dop853 does not have a 7th order embedded method (which is why it’s not in the name). DP8 is dop853.

Rational number

Optimality of stability vs leading truncation errors.

Ok, I think I understand that you’re using an implicit Runge-Kutta method that requires a nonlinear solve on 3 k variables at each time step for each variable in the ODE. I think.

So to answer my own questions,

What is the mathematical problem you’re trying to solve? Solving for coefficients of a implicit Runge-Kutta method on a nonlinear ODE.

What variable is nlsolve supposed to be searching over? Nsolve searches over k.

What is the dimension of that variable? The dimension of k is k=3n, where n is the dimension of the ODE.

What does the variable k represent? k represents Runge-Kutta coefficients for an implicit RK method on a nonlinear ODE. At any given time step the initial condition u is known and fixed and k is the solution to the implicit RK equations.

Is the search space u or k? The search space is k.

If all that is true, then the nonlinear solve should be within the Runge-Kutta timestepping function, inside the time-stepping loop, not at global scope.

But again I direct you to this bit of code.

function f!(fx, u)
   fx .= f(u, t, dt, k)
end

c1 = nlsolve(f!, k)

The input variable to your f! function is u. But when you call nlsolve you provide an initial guess for k, whose dimension is three times the dimension of u. This is likely the source of your dimension mismatch errors. To fix the errors you have to properly coordinate the search function and the search variables.

@ChrisRackauckas , I think @Vick wants to roll his own method in order to learn how to do things in Julia, and that it’s better to help him get his problem straight and his code working before delving into tricks for writing fast Julia code, determining the best integration algorithm, and solving the problem with the DifferentialEquations library. His problems with the Julia code are fundamental and not just suboptimal.

1 Like

I think you are right! I don’t understand how nlsolve is supposed to work. For me, I’m providing the function f which is outputting the correct system of non linear equations to solve and I’m also providing guesses in k. To me, providing the equations in the function of f and the guesses in terms of k, should have done it. Scipy fsolve does exactly that in Python.

I just don’t understand why we need to provide u to nlsolve, which is just the array of variables for any particular ODE function.

By the way, the code as I posted it, is not complete. It involved only the f function that had the As and Cs. The other piece of code is the complete Radau function that has the Bs and it is within this latter function that I’m calling nlsolve…

I just don’t understand why we need to provide u to nlsolve, which is just the array of variables for any particular ODE function

I think I misdirected you initially to provide u to nlsolve, because I misunderstood what the mathematical problem was. I think I understand correctly now that k is a set of implicit Runge-Kutta coefficients that have to be solved for at each time step, from a set of nonlinear equations involving f for fixed values of u and t So u and t are the parameters in the nonlinear solve and k is the unknown.

So the code snippet

function f!(fx, u)
   fx .= f(u, t, dt, k)
end

c1 = nlsolve(f!, k)

should be changed to

function f!(fk, k)
   fk .= f(u, t, dt, k)
end

solution = nlsolve(f!, k)
ksolution = solution.zero

And that code snippet should be inside the time-stepping loop, so that the u and t parameters can be set to their current values at any given time step.

2 Likes

That was it! Thanks a lot! My code with everything from scratch took 5 mins to run and solve for the van der Pol equation. But now having substituted nlsolve for my nrm functions (Newton and Jacobian matrix), my run is, as I was expecting; under 45 seconds.

Thanks again Mr. John Gibson.

2 Likes

Next you should look into Performance Tips · The Julia Language

Probably the main thing will be to preallocate vectors for return values and to remove temporary vectors that cause allocations within loops.

2 Likes