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

[Vick]

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:

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

[John_Gibson]

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

[Vick]

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:

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:

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?

[ChrisRackauckas]

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

[Vick]

I’d like to try to see for myself.

[John_Gibson]

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.

[Vick]

Have you tried your modifications?

Because it’s not working for me!

[John_Gibson]

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

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 …

[mkitti]

@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?

function f(u = [2, 0], t = 0, dt = 0.001, k = [1, 1, 1, 1, 1,1])
...

You probably need to rearrange the arguments.

[John_Gibson]

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

[Vick]

@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?

[Vick]

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]

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

[Vick]

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?

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

[zdenek_hurak]

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

[ChrisRackauckas]

@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:

SciML/OrdinaryDiffEq.jl/blob/master/src/perform_step/firk_perform_step.jl#L126-L143

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

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

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

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?

[ChrisRackauckas]

@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:

SciML/OrdinaryDiffEq.jl/blob/master/src/perform_step/firk_perform_step.jl#L126-L143

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.

[Vick]

@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++