Error for solving ODEs by an implicit RK methods when using NLsolve.jl

General Usage
,
1

Hi Guys!,

I’m trying to an ODE system using implicit RK method. Defining the following system
\dot q_i=p_i,\quad \dot p_i = -\frac{q_i}{||q||^3} and here q,p\in \mathbb R^3. I’m using this code

using NLsolve

function RK(f::Function,g::Function,t₀,tf,q₀,p₀,N)
    
 h = (tf-t₀)/N;
 t=t₀:h:tf; 

a = [0  0  0 ; 5/24  1/3  -1/24 ; 1/6  2/3  1/6];   
â = [1/6  -1/6  0 ; 1/6 1/3  0 ; 1/6  5/6  0];    

b = [1/6 , 2/3 , 1/6] ;    
    
q=[] ; push!(q,q₀);  p=[] ; push!(p,p₀); 

function k_stage!(F,ℓ,q₀,p₀)

F[1] = ℓ[1] .- g(q₀ .+ h * (â[1,1] * ℓ[4] + â[1,2] * ℓ[5] + â[1,3] * ℓ[6]));
F[2] = ℓ[2] .- g(q₀ .+ h * (â[2,1] * ℓ[4] + â[2,2] * ℓ[5] + â[2,3] * ℓ[6]));     
F[3] = ℓ[3] .- g(q₀ .+ h * (â[3,1] * ℓ[4] + â[3,2] * ℓ[5] + â[3,3] * ℓ[6]));
F[4] = ℓ[4] .- f(p₀ .+ h * (a[1,1] * ℓ[1] + a[1,2] * ℓ[2] + a[1,3] * ℓ[3]));
F[5] = ℓ[5] .- f(p₀ .+ h * (a[2,1] * ℓ[1] + a[2,2] * ℓ[2] + a[2,3] * ℓ[3]));        
F[6] = ℓ[6] .- f(p₀ .+ h * (a[3,1] * ℓ[1] + a[3,2] * ℓ[2] + a[3,3] * ℓ[3]));      
        
end

    
for i=1:N;
    
res = nlsolve((F,ℓ) -> k_stage!(F,ℓ,q₀,p₀), rand(6));

k = res.zero
    
q₀ = q₀ + (h/6) * (k[1]+2*k[2]+k[3]);
p₀ = p₀ + (h/6) * (k[4]+2*k[5]+k[6]); 

push!(q,q₀);  push!(p,p₀); 
        
end    
  return t, q, p; 
end

I opposite an error when applying the above code for the following problem

function f(p)   
return Qdot = [p[1] , p[2] , p[3]];
end

function g(q)
r = (√(q'*q))^3;
return Pdot = [q[1]/r , q[2]/r , q[3]/r];
end

t₀=0;  tf=100;   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]]; 

RK(f,g,t₀,tf,q₀,p₀,N)
2

I think the error is because your function g returns a 3-vector, which mean that ℓ[1] .- g(...) will also be a 3-vector, which you are subsequently trying to store in F[1], which will only store a scalar. I don’t follow what you’re trying to do well enough to suggest solutions though.

3

For a scalar function it works, let say for f(q)=q but didn’t for f(q)=[q[1]] even I have a scalar?