# Implicit numerical method using NLsolve.jl

**URL:** <https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525>\
**Category:** General Usage\
**Tags:** question, numerics, nlsolve, nonlinear, methods\
**Created:** [August 23, 2024, 10:34am UTC](https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525 "2024-08-23T10:34:00Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [August 23, 2024, 10:34am UTC](https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525/1 "2024-08-23T10:34:00Z")

</div>

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

```julia
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

```julia
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);

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 23, 2024, 10:59am UTC](https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525/2 "2024-08-23T10:59:50Z")

</div>

Note that Labatto tableaus have a lot of structure to them and it’s really not ideal that you solve them like this.

---

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [August 23, 2024, 11:06am UTC](https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525/3 "2024-08-23T11:06:05Z")

</div>

Thanks, but it’s an implicit and this is one way to do that. isn’t?

---

<div class="post-metadata">

**Author:** ![Andrew\_Swann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andrew_swann/32/32106_2.png) [@Andrew\_Swann](https://discourse.julialang.org/u/Andrew_Swann)\
**Post date:** [August 23, 2024, 11:13am UTC](https://discourse.julialang.org/t/implicit-numerical-method-using-nlsolve-jl/118525/4 "2024-08-23T11:13:33Z")

</div>

You get the same error if you just evaluate `g` on `ini_gauss`, so that needs to be corrected.
