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.
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:
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 ? Or if you really want a PCE symplectic method, the 16th order Gauss-Legendre is rather good and should outperform GJ quite handedly:
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.
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.
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.
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.
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.