Hi all. I am trying to simulate a geometric control which we proposed in our paper. Since it involves a differential equation evolving on a manifold. I thought of using manifold projection as mentioned in Callback Library · DifferentialEquations.jl . However, I am getting instability in solving the equation. Also not able to set up the NLSOLVEJL_SETUP(). I am including my code here. It would be really helpful if anyone could suggest what I could do differently here.

----------------------------------CODE--------------------------------------

using DifferentialEquations;

using OrdinaryDiffEq;

using Sundials;

using LinearAlgebra;

using StaticArrays;

using ODEInterface;

using ODEInterfaceDiffEq;

using LSODA;

function toFindLogRot!(RotMat)

```
trace=tr(RotMat);
println(trace)
arg = (trace-1)/2;
if arg >= 1
arg = 1
end
theta=acos(arg);
if theta != 0
omega=(1/(2*sin(theta)))*[RotMat[3,2]-RotMat[2,3];RotMat[1,3]-RotMat[3,1];RotMat[2,1]-RotMat[1,2]];
else
omega=[1;0;0];
end
return [omega' theta]'
```

end

----------------Manifold Part---------------

Since the rotation matrix has to satisfy transpose®*R =I, Can I write Like this??

function g(resid,R,p,t)

```
R = reshape(R,(3,3))
dR = transpose(R)*R - I
resid = reshape(dR,(1,:))
```

end

cb = ManifoldProjection(g,nlsolve=NLSOLVEJL_SETUP()) # How to setup nlsolve setup ??

function geometricControl!(dR,R,p,t)

```
R_if = reshape(R,(3,3))
p_if =[100*sin(t/20),100*sin(t/40),10];
p_if_dot = [5*cos(t/20);2.5*cos(t/40);0]
p_if_dotdot = [-0.25*sin(t/20);-(25/400)*sin(t/40);0]
p_fr = transpose(R_if)*(p_ir-p_if);
r = sqrt(p_fr[1]^2+p_fr[2]^2+p_fr[3]^2)
theta=acos( p_fr[3] / r);
phi=atan(p_fr[2] , p_fr[1]);
Rd = [ cos(phi)*cos(theta) cos(theta)*sin(phi) -sin(theta);
-sin(phi) cos(phi) 0;
cos(phi)*sin(theta) sin(phi)*sin(theta) cos(theta)];
R_if_dot = reshape(dR,(3,3))
p_fr_dot = transpose(R_if_dot)*(p_ir-p_if)-transpose(R_if)*p_if_dot
theta_dot = (-1/sqrt(p_fr[1]^2+p_fr[2]^2))*((r^2*p_fr_dot[3]-p_fr[3]*(p_fr[1]+p_fr[2]+p_fr[3]))/r^2)
phi_dot = (p_fr[1]*p_fr_dot[2]-p_fr[2]*p_fr_dot[1])/(sqrt(p_fr[1]^2+p_fr[2]^2))
Rd_dot = [-sin(phi)*cos(theta)*phi_dot-cos(phi)*sin(theta)*theta_dot cos(theta)*cos(phi)*phi_dot-
sin(theta)*sin(phi)*theta_dot -cos(theta)*theta_dot;
-cos(phi)*phi_dot -sin(phi)*phi_dot 0;
-sin(theta)*sin(phi)*phi_dot+cos(theta)*cos(phi)*theta_dot
sin(theta)*cos(phi)*phi_dot+cos(theta)*sin(phi)*theta_dot -sin(theta)*theta_dot]
er = toFindLogRot!(R_if*transpose(Rd))
eR = [0 -er[3] er[2];er[3] 0 -er[1];-er[2] er[1] 0];
omega = -beta * transpose(Rd)*er[4]*eR*Rd - transpose(Rd_dot)*Rd
R_if_dot = R_if*omega
dR .= reshape(R_if_dot,(1,:))
```

end

p=0

beta = 0.05

p_ir = [20,10,100];

tfinal = 2;

time_span = (0.0,tfinal);

R_if_0 = reshape(Matrix{Float64}(I,3,3),(1,:))

prob_di= ODEProblem(geometricControl!,R_if_0,time_span,p);

@time sol_di = solve(prob_di,CVODE_BDF(),callback=cb);

---------------------------ERROR----------------------

- Instability detected if nsolve is not used
- NLSOLVEJL_SETUP not defined

What I have tried —>

I have gone through this website Specifying (Non)Linear Solvers · DifferentialEquations.jl.

But code through an error CS and AD not defined (CS being chunk size, AD is Autodifferentiation (not sure about what AD mean)).

If possible please suggest what could be done about this.