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

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

I would do option 2.

12 个评论
显示 10更早的评论隐藏 10更早的评论

William Rose 2023-10-9

编辑：William Rose 2023-10-9

在 MATLAB Online 中打开

@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.

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.
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.

William Rose 2023-10-11

在 MATLAB Online 中打开

@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.

William Rose 2023-10-11

在 MATLAB Online 中打开

@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 2023-10-12

编辑：PONNADA 2023-10-12

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).

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