How to find largest Lyapunov exponent for the second order forced duffing oscillator by using the time series data [obtained by ODE45]?

Equation of the form:

dydt = [ x(2) ; f*sin(w*t)-c*x(2)-k1*x(1)-k3*x(1)^3]

where, [x(1), x(2)] =[ 0 , 0];

f = 0.5; w= 2 ; c = 0.01; k1 = 0.0213; k3 = 0.319

Answer 1

[William Rose]

@PONNADA,

If you have access to the Control Systems toolbox, do

lyapExp = lyapunovExponent(X,fs)

where X is the signal and Fs is its sampling rate.

Since the output from ode45 is not sampled at a constant rate, you will want to either

1. Pass a tspan argument to ode45() that includes the intermediate times at a constant rate, or
2. use interp1() to resample X to a constant rate, after ode45() finishes.

I would do option 2.

[William Rose]

@PONNADA,

Solve Duffing oscillator:

f = 0.5; w=2; c=0.01; k1=0.0213; k3=0.319;
tspan=[0,100];
y0=[0,0];
[t,y]=ode45(@(t,y) [y(2) ; f*sin(w*t)-c*y(2)-k1*y(1)-k3*y(1)^3], tspan, y0);

Resample at constant rate. I chose 10 Hz, after inspecting the uneven sampling intervals of the solution, the shortest of which was about 0.1 s.

Fs=10;                                      % sampling rate (Hz)
tInterp=tspan(1):1/Fs:tspan(2);             % time base for interpolation
yInterp=interp1(t,y(:,1),tInterp,'pchip');  % interpolate
figure; 
plot(t,y(:,1),'r.',tInterp,yInterp,'-b');   % plot signal
legend('raw','interp'); grid on; xlabel('Time (s)'); ylabel('y'); 

Estimate the Lyapunov exponent:

lyapExp = lyapunovExponent(yInterp,Fs)

lyapExp = 4.0693

Is the estimated L.E. different if you sample faster or slower, or if you choose a different initial condition, or if you increase or decrease the simulation duration by a factor of 2? What if you change the oscillator parameters to values that make the system periodic, instead of chaotic? Or make the system linear (k3=0)?

[PONNADA]

Dear William Rose, I am very glad for your kind response; I hope you are doing well. actually, In your given code, I need to make tspan =0:0.005:3500 to reach the steady state part of the solution very smoothly. Moreover, the solution for a given set of parameters is the periodic solution(period = 2*pi/w); in such a case, the Lyapunov exponent must be less than zero, but I am getting 40.369.

What should be the input for the Lyapunov exponent function?

Say, combinations: 1. LyapunovExponent( y(:,1), 1/0.005)

2. LyapunovExponent( y1_steady, 1/0.005) ; y1_steady = excluded transition part in y(:,1).

(or) neither 2, then give me further suggestions.

Thank you.

[William Rose]

@PONNADA, I plan to reply to your question later today, I hope.

[William Rose]

@PONNADA,

You wish to simulate for 3500 s and resample at 200 Hz. When I tried it in this window, I got an error: "Quitting, code takes longer than 55 s to run." or somehting like that. I think 200*3500=700,000 points is too many. So I tried it at 10 Hz, 20 Hz, and 50 Hz.

f = 0.5; w=2; c=0.01; k1=0.0213; k3=0.319;
tspan=[0,3500];     % Ponnada wants 3500 s 
y0=[0,0];           % initial condition
[t,y]=ode45(@(t,y) [y(2) ; f*sin(w*t)-c*y(2)-k1*y(1)-k3*y(1)^3], tspan, y0);

Plot the simulation results. Top panel shows all 3500 s. Bottom panel shows initial 50 s, middle 20 s, and final 20 s.

figure; 
subplot(211); plot(t,y(:,1),'-r');              % plot full length signal
grid on; xlabel('Time (s)'); ylabel('y'); 
subplot(234); plot(t,y(:,1),'-r');  
grid on; xlabel('Time (s)'); ylabel('y'); xlim([0 50])  % initial 50s
subplot(235); plot(t,y(:,1),'-r');   
grid on; xlabel('Time (s)'); xlim([1790 1810])  % middle 20 s
subplot(236); plot(t,y(:,1),'-r');   
grid on; xlabel('Time (s)'); xlim([3480 3500])  % final 20 s

If I resample at 200 Hz - the rate requetsted by @PONNADA. - I get an error message when I run the code in the online window, because the code takes too long to compute the Lyapunov exponent. Therefore I resample at 3 rates: 10 Hz, 20 Hz, and 50 Hz.

Fs10=10; Fs20=20; Fs50=50;          % sampling rates (Hz)
t10 =tspan(1):1/Fs10:tspan(2);      % 10 Hz time base
t20 =tspan(1):1/Fs20:tspan(2);      % 20 Hz time base
t50 =tspan(1):1/Fs50:tspan(2);      % 50 Hz time base
y10 =interp1(t,y,t10,'pchip'); % interpolate
y20 =interp1(t,y,t20,'pchip'); % interpolate
y50 =interp1(t,y,t50,'pchip'); % interpolate

Estimate the Lyapunov exponent from the full length signal, sampled at 10 Hz, 20 Hz, and 50 Hz, and display results on console:

lyapExp10 = lyapunovExponent(y10,Fs10);
lyapExp20 = lyapunovExponent(y20,Fs20);
lyapExp50 = lyapunovExponent(y50,Fs50); 
fprintf('Lyapunov exponent=%.2f at 10 Hz, %.2f at 20 Hz, %.2f at 50 Hz.\n',...
    lyapExp10,lyapExp20,lyapExp50);

Lyapunov exponent=1.55 at 10 Hz, 2.24 at 20 Hz, 5.98 at 50 Hz.

The results above are not what I expected, for several reasons.

1. The L.E. changes a lot when the sampling rate is changed. All the rates used are more than fast enough to capture the signal details, therefore I expected the LE would be about the same for all three sampling rates.
2. All three estimates of the LE are a lot more than 0. I expected LE to be approximatrely 0 for a periodic oscillation.

Let's try something even more simple: estimate the L.E. for a 1 Hz sine wave sampled at 10, 20, 50, 100 Hz. The signals x10, x20, etc. include the signal and its first derivative. I expect the LE to be approximately zero in all cases.

x10 =[cos(2*pi*[0:.10:10]'),-sin(2*pi*[0:.10:10]')];
x20 =[cos(2*pi*[0:.05:10]'),-sin(2*pi*[0:.05:10]')];
x50 =[cos(2*pi*[0:.02:10]'),-sin(2*pi*[0:.02:10]')];
x100=[cos(2*pi*[0:.01:10]'),-sin(2*pi*[0:.01:10]')];
LE10 =lyapunovExponent(x10,10);
LE20 =lyapunovExponent(x20,20);
LE50 =lyapunovExponent(x50,50);
LE100=lyapunovExponent(x100,100);
fprintf('LE10=%.2f, LE20=%.2f, LE50=%.2f, LE100=%.2f.\n',LE10,LE20,LE50,LE100)

LE10=0.56, LE20=-0.14, LE50=0.23, LE100=5.27.

Interesting. The first three estimates are not bad, but the last one is quite unexpected. Perhaps function lyapunovExponent() works poorly if the signal is very much oversampled.

[PONNADA]

Dear William, your care and patience are evident in your reply. Your unrewarded work always gives you the best. Moreover, thank so much for your kind reply.

[William Rose]

@PONNADA, thank you for your kind words.

I am not entirely certain whether I should pass y1(t) and y2(t)=dy1/dt to lyapunovExponent(), or pass only y1(t). The function gives answers with either choice. In the Matlab help for lyapunovExponent(), they show passing only one component of the phase space for the 3-dimensional Lorenz attractor.

The method used in lyapunovExponent() is described birefly in the Matlab help. The original description of the algorithm is attached (Rosenstein et al., 1993). The senior author is Carlo De Luca. I am familiar with some of his work.

I think 700,000 points (3500 seconds times 200 samples per second) is far more than necessary. The examples in the original paper use 500 to 8000 points. The algorithm may not work well if the sampling interval is too short. This explains the results we a-observed above, when we sampled a 1 Hz sine wave at 100 Hz, and got unreliable estimates of the Lyapunov exponent. SInce the Duffing oscillator with the parameters chosen has a final frequency of about 0.3 Hz, a 5 Hz sampling rate gives 15-20 samples per cycle. Rosenstein et al. recommend n to 10n points per orbit of the attractor, where n=dimension of the system; n=2 for our Duffing oscillator. Therefore we want 2 to 20 points per orbit, and we have that by sampling at 5 Hz.

Simulate Duffing oscillator:

f = 0.5; w=2; c=0.01; k1=0.0213; k3=0.319;
tspan=[0,100];     % Ponnada wants 3500 s 
y0=[0,0];           % initial condition
[t,y]=ode45(@(t,y) [y(2) ; f*sin(w*t)-c*y(2)-k1*y(1)-k3*y(1)^3], tspan, y0);

Resample at 5 Hz:

Fs=5;                           % sampling rates (Hz)
t5 =tspan(1):1/Fs:tspan(2);     % 5 Hz time base
y =interp1(t,y,t5,'pchip');     % interpolate

Estimate lag

dim=2;
[~,lag] = phaseSpaceReconstruction(y(:,1),[],dim)

lag = 7

Estimate L.E. at many expansion steps:

eRange = 100;
lyapunovExponent(y,Fs,lag,dim,'ExpansionRange',eRange)

One is supposed to use the movable selector lines on the plot above to select a range of epansion steps - see Help. But I cannot run the script in Matlab since I lack the toolbox, and the plot above is not adjustable. The plot above also does not look like the exa[les in the Help, so I am not sure how I would choose the expension step range. I will try the range 40 to 90, since that looks relatively linear.

Kmin=40; Kmax=90;

Estimate the Lyapunov exponent from the full length signal, sampled at 5 Hz, and display results on console:

lyapExp5 = lyapunovExponent(y,Fs,lag,dim,'ExpansionRange',[Kmin,Kmax]);
fprintf('Lyapunov exponent=%.3f at 5 Hz.\n',lyapExp5);

Lyapunov exponent=0.054 at 5 Hz.

I would try to use lypaunovExponent() to reproduce results forom Rosenstein et al. 1993, and from other sources. This would help improve your ability to use the function effectively.

[PONNADA]

Dr.William, I believes, I can get output from you. Thank you.

[William Rose]

@PONNADA,

Rosenstein et al. 1993, in the article whicich attached to my previous comment, say on p. 122:

"For each system, the initial point was chosen near the attractor and the transient points were discarded. In all cases, the x-coordinate time series was used to reconstruct the dynamics."

To me, the statement "the x-coordinate time series was used" answers the question about whether to pass the two-column array y or just the vector y(:.1): I should passonly the vector y(:,1). Also, I had not been chopping off the initial transient, and Rosenstein says I should have chopped off the initial transient.

Therefore I will simulate out to t=3500 s, as you initially recommended, and I will discard the data from t=0 to 1500.

Simulate Duffing oscillator:

f = 0.5; w=2; c=0.01; k1=0.0213; k3=0.319;
tspan=[0,3500];     % Ponnada wants 3500 s 
y0=[0,0];           % initial condition
[t,y]=ode45(@(t,y) [y(2) ; f*sin(w*t)-c*y(2)-k1*y(1)-k3*y(1)^3], tspan, y0);

Resample y(:,1) at 5 Hz, starting at time 1500 s:

Fs=5;                           % sampling rate (Hz)
t5=1500:1/Fs:tspan(2);          % 5 Hz time base
y5=interp1(t,y,t5,'pchip');     % interpolate

Plot results: y1 from 1500 to 3500 s. Plot y1 versus time and phase plots for short segments.

subplot(311); plot(t5,y5(:,1)); grid on; xlabel('Time (s)'); ylabel('y_1') 
subplot(334); plot(t5,y5(:,1)); grid on; xlabel('Time'); ylabel('y_1'); xlim([1500,1540])
subplot(335); plot(t5,y5(:,1)); grid on; xlabel('Time'); xlim([2500,2540])
subplot(336); plot(t5,y5(:,1)); grid on; xlabel('Time'); xlim([3460,3500])
subplot(337); plot(y5(t5<1540,1),y5(t5<1540,2)); 
grid on; ylabel('dy_1/dt'); xlabel('y_1')
subplot(338); plot(y5(t5<2540,1),y5(t5<2540,2)); 
grid on; xlabel('y_1')
subplot(339); plot(y5(t5>3460,1),y5(t5>3460,2)); 
grid on; xlabel('y_1')

The plots above show that the system has acheived a steady state oscillaiton by the time t>=1500.

Estimate lag:

dim=2;
[~,lag] = phaseSpaceReconstruction(y5(:,1),[],dim)

lag = 4

Estimate L.E. at many expansion steps:

eRange = 100;
lyapunovExponent(y5(:,1),Fs,lag,dim,'ExpansionRange',eRange)

One is supposed to use the movable selector lines on the plot above to select a range of expansion steps - see Help. The plot above also does not look like the examples in the Help, nor does it look like the examples in Rosenstein et al, 1993. Therefore I am not sure how I would choose the expansion step range. I will try the range 5 to 50, but really, I expect the slope of the line in the plot above (which is the estimate of the L.E.) is about 0, for any range of Kmin,Kmax, from 1 to 100..

Kmin=5; Kmax=50;

Estimate the Lyapunov exponent from the full length signal, sampled at 5 Hz, and display results on console:

lyapExp5 = lyapunovExponent(y5(:,1),Fs,lag,dim,'ExpansionRange',[Kmin,Kmax]);
fprintf('Lyapunov exponent=%.3f at 5 Hz.\n',lyapExp5);

Lyapunov exponent=0.005 at 5 Hz.

I tried it with Fs=10 Hz and the results were similar.

[PONNADA]

Dr.William, I am delighted to receive your kind and immediate response. It makes me more interested and enjoyable to learn about LapunovExponent. I gave the reference Fig.1 below[Ref:doi:https://doi.org/10.1016/j.chaos.2023.113820 , Fig .14], in this subplot(121), ϕ Vs λmeans w Vs y1_L(1500:3500, 1)[ y1_L(1500:3500,1) = local maxima and minima in t Vs y(1500:3500,1) ] in my case( and also subplot 336 in your latest comment). So, if I convert subplot(336) given in the older comment to w Vs local maxima and minima in y(1500:3000,1), I got Fig. 2. Hence, this solution is a periodic means Lyapunov exponent must be less than zero as referred in Fig. 01. But why it is >0 (0.000356563)? Please help me in this regard( to find the regions of chaos and periodic solutions in the particular range of 'w', using the Lyapunov exponent).

%---------------------------------------------------------------------------------------------------------------------------------------------

Fig.1.

%---------------------------------------------------------------------------------------------------------------------------------------------

Fig. 2.

%---------------------------------------------------------------------------------------------------------------------------------------------

[William Rose]

@PONNADA,

You asked "why it is >0 (0.000356563)?"

I think the esitmate of the LE is >0 because the method of estimation is not perfect. Remember that the method of esitmation used is to compute the slope of the bestf-it line thorugh the pointsin this plot (which I have copied from my comment above):

The user specifies the range of expansion steps (=x-values) to use when computing the slope. If you do not specify a range, the full range is used, which is K=1..100 in the example above. The slope was L.E.=0.00035 for this range. When I selected the range K=5:50, I got L.E.=0.005. The article by Rosenstein et al., 1993, does not give clear guidance on how to select the range of expansion steps.

[PONNADA]

Thank you Dr. William for your valid information and having a nice cummunication. I rejoice in your goodness.

[William Rose]

👍. Thank you for your kind comments, @PONNADA.