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?