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, p * (1 - u^2) * u - u] 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?