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

**URL:** <https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349>\
**Category:** New to Julia\
**Tags:** package\
**Created:** [July 17, 2022, 8:23am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349 "2022-07-17T08:23:31Z")\
**Posts on this page:** 13\
**Page:** 2

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 18, 2022, 12:14pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/21 "2022-07-18T12:14:20Z")

</div>

> [@ChrisRackauckas](#):
>
> 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.

By the way, I don’t see Dopri 8(7) with 13 stages listed in differential Equations.jl Why is that? [article](http://peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK8/RKcoeff8d_1.pdf)

There is another algorithm not listed also, but which I really needed. The Gauss-Jackson 8th Order for orbit propagation - [article](https://www.researchgate.net/publication/333300756_Accuracy_and_Efficiency_Comparison_of_Six_Numerical_Integrators_for_Propagating_Perturbed_Orbits)  
and [article2](https://drum.lib.umd.edu/bitstream/handle/1903/2202/2004-berry-healy-jas.pdf;sequence=7)

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 18, 2022, 1:55pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/22 "2022-07-18T13:55:39Z")

</div>

> [@Vick](#):
>
> 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.

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](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757) and follow its recommendations.

---

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 18, 2022, 3:03pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/23 "2022-07-18T15:03:13Z")

</div>

Ok  
Here is my code for nlsolve:

```julia
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]

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 18, 2022, 10:07pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/24 "2022-07-18T22:07:16Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 19, 2022, 2:58am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/25 "2022-07-19T02:58:44Z")

</div>

> [@Vick](#):
>
> By the way, I don’t see Dopri 8(7) with 13 stages listed in differential Equations.jl Why is that? [article](http://peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK8/RKcoeff8d_1.pdf)

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:

> <https://github.com/SciML/DiffEqDevTools.jl/blob/master/src/ode_tableaus.jl>

and documented here:

[https://docs.sciml.ai/dev/modules/DiffEqDevDocs/internals/tableaus/](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](https://arxiv.org/abs/2207.08135)

 ![Screenshot 2022-07-18 221818](https://global.discourse-cdn.com/julialang/original/3X/1/e/1e7086009a04acda42e3ea39e70344d693df5d56.png)

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.

> [@Vick](#):
>
> There is another algorithm not listed also, but which I really needed. The Gauss-Jackson 8th Order for orbit propagation - [article](https://www.researchgate.net/publication/333300756_Accuracy_and_Efficiency_Comparison_of_Six_Numerical_Integrators_for_Propagating_Perturbed_Orbits)  
> and [article2](https://drum.lib.umd.edu/bitstream/handle/1903/2202/2004-berry-healy-jas.pdf;sequence=7)

Ehh, it’s in the todo list:

> <https://github.com/SciML/OrdinaryDiffEq.jl/issues/331>
>
> Gauss-Jackson: https://drum.lib.umd.edu/bitstream/handle/1903/2202/2004-berry-he…aly-jas.pdf?sequence=7
> https://www.researchgate.net/publication/333300756\_Accuracy\_and\_Efficiency\_Comparison\_of\_Six\_Numerical\_Integrators\_for\_Propagating\_Perturbed\_Orbits
> 
> Berry-Healy: https://pdfs.semanticscholar.org/2166/2dcd9b7ece9f4abbb2e52154ae7a840ed0d6.pdf or https://pdfs.semanticscholar.org/4701/d1c3c9b1101c3e23b4472d4895aa5d7f9a51.pdf
> 
> Stormer-Cowell: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.573.319&rep=rep1&type=pdf and https://www.sciencedirect.com/science/article/pii/S037704270400250X and https://www.sciencedirect.com/science/article/pii/0377042794001020
> 
> Hairer I has sections on the implicit and explicit Stromer methods as well.
> 
> Interesting related result: https://researchspace.auckland.ac.nz/bitstream/handle/2292/4980/449.pdf?sequence=1

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 [https://diffeq.sciml.ai/stable/solvers/dynamical\_solve/#Symplectic-Integrators](https://diffeq.sciml.ai/stable/solvers/dynamical_solve/#Symplectic-Integrators) ? Or if you really want a PCE symplectic method, the 16th order Gauss-Legendre is rather good and should outperform GJ quite handedly:

> **[GitHub - SciML/IRKGaussLegendre.jl: Implicit Runge-Kutta Gauss-Legendre 16th...](https://github.com/SciML/IRKGaussLegendre.jl)**
>
> Implicit Runge-Kutta Gauss-Legendre 16th order (Julia) - GitHub - SciML/IRKGaussLegendre.jl: Implicit Runge-Kutta Gauss-Legendre 16th order (Julia)

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.

---

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 19, 2022, 10:49am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/26 "2022-07-19T10:49:43Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 19, 2022, 11:07am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/27 "2022-07-19T11:07:35Z")

</div>

> [@ChrisRackauckas](#):
>
> `constructDormandPrice8()` is the method you refer to. But in no benchmark does that method do well.

The coefficients you have here is different from that of Hairer: [link](http://www.unige.ch/~hairer/prog/nonstiff/dop853.f)

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](https://en.wikipedia.org/wiki/Leapfrog_integration#Yoshida_algorithms) and Dopri 8(7) 13 for nbody problem.

I am interested in your Yoshida 6th order.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 19, 2022, 12:20pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/28 "2022-07-19T12:20:44Z")

</div>

> [@Vick](#):
>
> The coefficients you have here is different from that of Hairer: [link](http://www.unige.ch/~hairer/prog/nonstiff/dop853.f)

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`.

> <https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/tableaus/high_order_rk_tableaus.jl#L667-L758>

> [@Vick](#):
>
> What does “//” mean in Julia?

Rational number

> [@Vick](#):
>
> I do have Efficient Verner 9th Order. What’s the difference between the “Robust” and “Efficient” one?

Optimality of stability vs leading truncation errors.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 19, 2022, 1:08pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/29 "2022-07-19T13:08:51Z")

</div>

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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 19, 2022, 2:16pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/30 "2022-07-19T14:16:49Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 19, 2022, 2:40pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/31 "2022-07-19T14:40:06Z")

</div>

> 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

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

c1 = nlsolve(f!, k)

```

should be changed to

```julia
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.

---

<div class="post-metadata">

**Author:** ![Vick](https://avatars.discourse-cdn.com/v4/letter/v/aeb1de/32.png) [@Vick](https://discourse.julialang.org/u/Vick)\
**Post date:** [July 19, 2022, 3:52pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/32 "2022-07-19T15:52:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 20, 2022, 3:35pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/33 "2022-07-20T15:35:06Z")

</div>

Next you should look into [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/)

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

[Previous page](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349.md?page=1)