I would like to implement Lobatto IIIA-IIIB which is an implicit method. To do so, a non-linear equations need to be solved. Then I used NLsolve.jl for my problem
using NLsolve
function NuMethod(fV::Function,fT::Function,x₀,y₀,h,N)
x=[] ; push!(x,x₀); y=[] ; push!(y,y₀);
A = [0 0 0 ; 5/24 1/3 -1/24 ; 1/6 2/3 1/6];
B = [1/6 -1/6 0 ; 1/6 1/3 0 ; 1/6 5/6 0];
b = [1/6 , 2/3 , 1/6] ;
g(K,L) = vcat(
K[1] - fV(y₀+h*(B[1,1]*L[1]+B[1,2]*L[2]+B[1,3]*L[3])),
K[2] - fV(y₀+h*(B[2,1]*L[1]+B[2,2]*L[2]+B[2,3]*L[3])),
K[3] - fV(y₀+h*(B[3,1]*L[1]+B[3,2]*L[2]+B[3,3]*L[3])),
L[1] - fT(x₀+h*(A[1,1]*K[1]+A[1,2]*K[2]+A[1,3]*K[3])),
L[2] - fT(x₀+h*(A[2,1]*K[1]+A[2,2]*K[2]+A[2,3]*K[3])),
L[3] - fT(x₀+h*(A[3,1]*K[1]+A[3,2]*K[2]+A[3,3]*K[3])))
g(v) = g(v[1:3],v[4:6]); # a function of a single argument
for i in 1:N;
ini_gauss = vcat(fV(y₀),fV(y₀),fV(y₀),fT(x₀),fT(x₀),fT(x₀))
r = nlsolve(g, ini_gauss)
K = r.zero
x₀ = x₀ + h * ( b[1] * K[1] + b[2] * K[2] + b[3] * K[3]) ;
y₀ = y₀ + h * ( b[1] * k[4] + b[2] * K[5] + b[3] * K[6]) ;
push!(x,x₀) ; push!(y,y₀); # stock my solutions
end
return stack(x, dims=1), stack(y, dims=1); # reformulate solutions to vectors
end
I would really like to use this code for any functions fV and fT single or multidimensional case. However, I always opposite the error that MethodError: no method matching +(::Vector{Float64}, ::Float64) For element-wise addition, use broadcasting with dot syntax: array .+ scalar
This is my test probelm
function fKeplerT(p) # fKeplerV = vector field to potentioal energy
return Qdot = [p[1] , p[2] , p[3]];
end
function fKeplerV(q)
r = (√(q'*q))^3;
return Pdot = [-q[1]/r , -q[2]/r , -q[3]/r];
end
h=0.01; N=1000; x₀ = [0.8 , 0.6 , 0 , 0 , 1 , 0.5];
q₀=[x₀[1], x₀[2], x₀[3]]; p₀=[x₀[4], x₀[5], x₀[6]];
q,p = Lobatto_IIIA_IIIB(fKeplerT,fKeplerV,q₀,p₀,h,N);