# 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:** 20
**Page:** 1
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [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/1 "2022-07-17T08:23:31Z")
Hi,
I’m new to Julia. I have successfully coded a runge kutta implicit method in Julia on Windows using VS code.
I have built everything from scratch. However, the running time is about 5 minutes. So I want to use NLsolve to find root faster for the implicit solver
My system of ODEs to solve within the framework of the implicit solver are as follows:
```julia
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
k = ones(3*m)
c1 = nlsolve(f, k)
```
However I’m receiving errors:
> ERROR: MethodError: no method matching f(::Vector{Float64}, ::Vector{Float64})
> Closest candidates are:
> f(::Any, ::Any, ::Any, ::Any) at c:\Users\Vickram\Downloads\juli\rad5.jl:6
> f(::Any, ::Any, ::Any, ::Any, ::Any) at
t is the time (so for one instance of the solver, the t = 0); (t = (0, 3000) )
dt is the timestep (in this example: 0.001)
u is the initial conditions (in this example: [2, 0] )
The nlsolve should solve for one instance only as such:
t = 0, dt = 0.001, u = [2, 0]; and k = [1, 1, 1, 1, 1,1]
Please help me out. Thanks
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [July 17, 2022, 3:15pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/2 "2022-07-17T15:15:35Z")
Welcome to the Julia community! I suspect the fundamental problem here is that your `f` function does not fit the function signature that `nlsolve` expects. To find a zero of the mathematical function f(x), the Julia function should take the form `f!(fx, x)`, where `x` is the input vector and `fx` is a mutating storage vector for the output value f(x). For example,
```julia
julia> using NLsolve
julia> function f!(fx, x)
fx[1] = (x[1]+3)*(x[2]^3-7)+18
fx[2] = sin(x[2]*exp(x[1])-1)
end
f! (generic function with 1 method)
julia> solution = nlsolve(f!, [0.1; 0.2])
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [0.1, 0.2]
* Zero: [3.771277437827691e-13, 1.0000000000009224]
* Inf-norm of residuals: 0.000000
* Iterations: 5
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 6
* Jacobian Calls (df/dx): 6
julia> x = solution.zero
2-element Vector{Float64}:
3.771277437827691e-13
1.0000000000009224
```
To adapt your runge-kutta integrator to work this way, you need to write an `f!(fx,x)` function like that that calls your `f` with fixed `dt,k` parameters and sets `fx` to the output of your `f`, whatever that is. I tried to run your code and figure out what it’s doing, but it fails due to `m` being undefined.
I hope that the above suggestion is enough to get you on the right path.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 17, 2022, 4:02pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/3 "2022-07-17T16:02:28Z")
Hi, thanks for replying! The m is the same as the n in the f function. The m is actually in another function. I extracted the parts with m and k =ones…The k should have 6 ones in vector form: 3 ones for each variable in the original function:
```julia
function fu(u, t)
p = 1000
du = [
u[2],
p * (1 - u[1]^2) * u[2] - u[1]]
return du
end
```
This is the van der Pol Equation.
so:
```julia
m = length(u)
t = 0.0 # begin
dt = 0.001 # timestep
tn = 3000 # end
u = [2.0, 0.0]
```
Do you still think I need to modify my code?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 17, 2022, 4:45pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/4 "2022-07-17T16:45:30Z")
> [@Vick](#):
>
> I have built everything from scratch. However, the running time is about 5 minutes. So I want to use NLsolve to find root faster for the implicit solver
Note that the algorithm you’re using here is not efficient, so you could optimize the code a bit but it will still be slow.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 17, 2022, 4:49pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/5 "2022-07-17T16:49:36Z")
I’d like to try to see for myself.
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [July 17, 2022, 4:57pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/6 "2022-07-17T16:57:38Z")
Yes, you need to modify your code. Either change your `f` to match the function signature that NLsolve expects, or write another function with the right signature that calls your `f` with some fixed integration parameters.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 17, 2022, 6:41pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/7 "2022-07-17T18:41:29Z")
Have you tried your modifications?
Because it’s not working for me!
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [July 17, 2022, 6:49pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/8 "2022-07-17T18:49:09Z")
No, I haven’t. It’s too much work for me to figure out what you mean without a mathematical specification of the problem and code that doesn’t run. Maybe I can give you a complete working version of something similar for reference.
> [@Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …
---
### Author:  [@mkitti](https://discourse.julialang.org/u/mkitti)
#### Post date: [July 17, 2022, 7:00pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/9 "2022-07-17T19:00:55Z")
> [@Vick](#):
>
> t = 0, dt = 0.001, u = [2, 0]; and k = [1, 1, 1, 1, 1,1]
How about you start with something lkke this?
```julia
function f(u = [2, 0], t = 0, dt = 0.001, k = [1, 1, 1, 1, 1,1])
...
```
You probably need to rearrange the arguments.
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [July 17, 2022, 7:05pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/10 "2022-07-17T19:05:22Z")
Here’s a complete working example. My `f` is like your ODE integrator, with a couple parameters `a,b`. I then define an `f!(fx, x)` function that calls `f` with fixed parameters and has the right signature for `nlsolve`.
```julia
julia> function f(x, a, b)
[(x[1]+3)*(x[2]^3-7) + a; sin(x[2]*exp(x[1]) - b)]
end
f (generic function with 1 method)
julia> function f!(fx, x)
fx .= f(x, 18, 1)
end
f! (generic function with 1 method)
julia> using NLsolve
julia> solution = nlsolve(f!, [0.1; 0.2])
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [0.1, 0.2]
* Zero: [3.771277437827691e-13, 1.0000000000009224]
* Inf-norm of residuals: 0.000000
* Iterations: 5
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 6
* Jacobian Calls (df/dx): 6
julia> x = solution.zero
2-element Vector{Float64}:
3.771277437827691e-13
1.0000000000009224
```
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 5:59am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/11 "2022-07-18T05:59:06Z")
> [@John\_Gibson](#):
>
> `solution.zero`
In nlsolve ( … ) You separated the numbers with a semi-colon; what is the difference between a semi-colon and a comma in Julia?
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 6:12am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/12 "2022-07-18T06:12:00Z")
I cannot do that! 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) …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]
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 6:15am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/13 "2022-07-18T06:15:13Z")
I don’t think so: the same code on Python took about 1 hour (we are talking here about the Van der Pol equation) with everything built from scratch except for the fsolve from scipy.
My code in Julia is virtually the same, except that I’m not using any third party solver for nonlinear equations. I’m using my own code for Newton’s method. And this at present is taking about 5 minutes give or take 1 minute.
So scraping my own code for optimize solving, and using NLsolve, I’m expecting this to run under 1 minute.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 6:31am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/14 "2022-07-18T06:31:07Z")
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](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Radau_IIA_methods)
Therefore I have set up my code as follows:
ODE function for van der Pol equation:
```julia
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):
```julia
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:
```julia
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):
```julia
m = length(u)
k = ones(3*m)
```
And now the solver:
```julia
#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?
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 6:37am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/15 "2022-07-18T06:37:10Z")
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.
---
### Author:  [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)
#### Post date: [July 18, 2022, 8:42am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/16 "2022-07-18T08:42:34Z")
[https://docs.julialang.org/en/v1/base/punctuation/](https://docs.julialang.org/en/v1/base/punctuation/)
[https://juliateachingctu.github.io/Julia-for-Optimization-and-Learning/dev/lecture\_02/arrays/](https://juliateachingctu.github.io/Julia-for-Optimization-and-Learning/dev/lecture_02/arrays/)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 18, 2022, 10:31am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/17 "2022-07-18T10:31:47Z")
> [@Vick](#):
>
> I don’t think so: the same code on Python took about 1 hour (we are talking here about the Van der Pol equation) with everything built from scratch except for the fsolve from scipy.
>
> My code in Julia is virtually the same, except that I’m not using any third party solver for nonlinear equations. I’m using my own code for Newton’s method. And this at present is taking about 5 minutes give or take 1 minute.
>
> So scraping my own code for optimize solving, and using NLsolve, I’m expecting this to run under 1 minute.
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.
> **[coimbra.pdf](https://www.unige.ch/~hairer/preprints/coimbra.pdf)**
>
> 1233.85 KB
So actual implicit RK code would never use the algorithm you’re using here.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 11:16am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/18 "2022-07-18T11:16:10Z")
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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 18, 2022, 11:23am UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/19 "2022-07-18T11:23:18Z")
> [@Vick](#):
>
> 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.
The usually implementation, Hairer’s `radau` and `radau5` implementation in Fortran, both use this trick.
> [@Vick](#):
>
> 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?
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)](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](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.
---
### Author:  [@Vick](https://discourse.julialang.org/u/Vick)
#### Post date: [July 18, 2022, 12:01pm UTC](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349/20 "2022-07-18T12:01:17Z")
> [@ChrisRackauckas](#):
>
> Yes, in DifferentialEquations.jl. The code I linked is the Julia code that uses these tricks:
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++
[Next page](https://discourse.julialang.org/t/nlsolve-for-find-root-for-a-system-of-ode-for-an-implicit-runge-kutta-method/84349.md?page=2)