# Using spectral methods and optimization to find period of a set of nonlinear autonomous differential equations

hi all, I’m very new to Julia. I’m following these lectures to solve optimization problem.
I will be providing a toy problem (2 DoF) which mimics my actual problem (12 DoF).
I have a set of nonlinear non-autonomous differential equations . (dx, dy). For some value of parameter “a” this system has a periodic solution, with a period “T”.

‘’'julia

``````  function dzdt(x,y,t,a)
u    =  (1-(x^2)/4)-y^2
dx   =  -4*y+x*u  + a*x*exp((10^-3)*t)
dy   =   x+y*u    + a*y*exp((10^-2)*t)
[ dx, dy ]
end
``````

‘’’
There might be many methods to solve these kind of problems. I’m using collocation and spectral methods for periodic system to find periodic solution of given system.
Key idea

1. let “x(t), y(t)” be periodic solution to system “dzdt” for some “a”. We can easy represent a some random periodic function in Fourier series. Instead of continuous time here we deal with discrete time.
1.a. If we can represent solutions-> ( x(t), y(t) ) in terms of some Fourier series → ( F1(t), F2(t) ).
1.b. We can take time derivative of these solutions → ( d (x(t), y(t) ) /dt ) which should be our given set of differential equations.
1.c. Difference between guess periodic solution (actually the derivative of the guess solution) and actual differential equations should be zero. Then we can say that we found the periodic solution.

2. We construct this solution by assuming for some discrete time points (t = t1,t2,t3…tn) and at these time points ( x(t) , y(t) ) pass through some collocation points ( ( x1,y1), (x2,y2), (x3,y3),… (xn,yn) ). These (xi’s, yi’s) are our guess points.

3. We have matrix called "D Matrix for Fourier " . when this “D” matrix is pre-multiplied with the vector (v) it gives the approximate solution of the derivative of the (v) vector.
Here is the D Matrix

‘’'julia

``````function D(N::Integer)
h=2*pi/N;
ns=range(1,N-1,step=1);
col1(nps)=0.5*((-1)^nps)/sin(nps*h/2);
col=[0,col1.(ns)...];
row=[0,col[end:-1:2]...];
Toeplitz(col,row)
end
``````

‘’’
Where “N” is the number of discrete time points.

Discrete points in time. Since I don’t know the period of system I will take some high value as final time “tfinal” . Here I will construct some initial guess solution points for x(t),y(t).

‘’’ julia

``````tfinal     =    1.1*pi;
tpoints    =  collect( range(1,N,step=1) )*tfinal/N;
xguess     =  sin.( (2*pi/tfinal)*tpoints  )*2.0;
yguess     =  -sin.( (2*pi/tfinal)*tpoints )*0.5;
``````

‘’’
Here I’m writing a function to find the difference of assumed derivative solution and actual system.
I have written this function such that this can be used as a constraint but not as a objective function.

My objective function could be “0” . as of now.

‘’’ julia

``````  function dxlist(xs,ys,tf,a)
ns    =  2;
ts        =  collect( range(1,N,step=1) )*tf/N;
xytsZip   =  zip(xs,ys,ts);
dxD0      = [ dzdt(x,y,t,a) for (x,y,t) in xytsZip ];
dxD       = reduce(hcat, dxD0)';
xlyl      = reshape([xs;ys],N,ns);
dxF       = (Dmat*xlyl)*(2.0*pi/tf);
err       = dxD - dxF;

[ vcat(err'...)  .-10^(-10)  ;  -vcat(err'...)  .+10^(-10) ]
end

function cons(x)
# nps=(length(x[:])-1) ÷ 2;
tf=x[end-1];
a=x[end];
xs1=x[1:N];
ys1=x[N+1:2*N];
dxlist(xs1,ys1,tf,a)

end
``````

‘’’
initial guess points as a vector
‘’’ julia

``````a0   =  10^-3;

x0  =  vcat( [ xguess; yguess; [ tfinal , a0] ] );

xlower1=push!( -3.0*ones(2*N), pi );
xlower=push!(xlower1,-10^-3);
xupper1=push!(3.0*ones(2*N),1.5*pi);
xupper=push!(xupper,10^-3);
consLower=-ones(4*N)*Inf;
consUpper=zeros(4*N);
``````

‘’’
Here im calling optimization functions

‘’'julia

``````obj(x)=0.0;
model=ADNLPModel(obj, x0, xlower, xupper, cons, consLower, consUpper)
output=ipopt(model)
xstar=output.solution
fstar=output.objective
``````

‘’’
Packages used are

‘’'julia

``````  using LinearAlgebra
using DiffEqOperators, ForwardDiff
using ApproxFun
using FFTW
using ToeplitzMatrices
using ModelingToolkit
using DifferentialEquations
using NLPModels
using NLPModelsIpopt
using DCISolver
using JSOSolvers
``````

‘’’

Whole code in one place.
‘’’ julia

``````using LinearAlgebra
using DiffEqOperators, ForwardDiff
using ApproxFun
using FFTW
using ToeplitzMatrices
using ModelingToolkit
using DifferentialEquations
using NLPModels
using NLPModelsIpopt
using DCISolver
using JSOSolvers

# Number of collocation points
N=3

function Dmatrix(N::Integer)
h=2*pi/N;
ns=range(1,N-1,step=1);
col1(nps)=0.5*((-1)^nps)/sin(nps*h/2);
col=[0,col1.(ns)...];
row=[0,col[end:-1:2]...];
D=Toeplitz(col,row)
end

Dmat=Dmatrix(N);

function dzdt(x,y,t,a)
u=(1-(x^2)/4)-y^2;
dx=-4*y+x*u+a*x;
dy=x+y*u+a*y;
[dx,dy]
end
# initial guess
tfinal=1.1*pi;
tpoints=collect(range(1,N,step=1))*tfinal/N;
xguess=sin.((2*pi/tfinal)*tpoints)*2.0
yguess=-sin.((2*pi/tfinal)*tpoints)*0.5

function dxlist(xs,ys,tf,a)
nstates=2
ts=collect(range(1,N,step=1))*tf/N;
xytsZip=zip(xs,ys,ts);
dxD0=[dzdt(x,y,t,a) for (x,y,t) in xytsZip];
dxD=reduce(hcat, dxD0)';
xlyl=reshape([xs;ys],N,nstates);
dxF=(Dmat*xlyl)*(2.0*pi/tf);
err=dxD-dxF;
[vcat(err'...).-10^(-10);-vcat(err'...).+10^(-10)]
end

function cons(x)
tf=x[end-1];
a=x[end];
xs1=x[1:N];
ys1=x[N+1:2*N];
dxlist(xs1,ys1,tf,a)

end
a0=10^-3;
x0=vcat([xguess;yguess;[tfinal,a0]]);

obj(x)=0.0;

xlower1=push!(-3*ones(2*N),pi);
xlower=push!(xlower1,-10^-3)
xupper1=push!(3*ones(2*N),1.5*pi);
xupper=push!(xupper,10^-3)
consLower=-ones(4*N)*Inf;
consUpper=zeros(4*N)

# println("constraints vector = ",cons(x0))