0
$\begingroup$

For now, I have a coupled system with 5 variables and use the Runge-Kutta method to integrate.

sigma=9.95;k=28;b=8/3;C1=0.1;C2=1;Od=1;Om=10;Sm=10;
Ss=1;Spd=10;Sigma=100;C3=0.01;C4=0.01;C5=1;C6=0.001;

ntime=1000;
X=zeros(5,ntime);
dt=0.001; T=dt*ntime;
X(:,1)=[0,1,0,0,0]';
for t = 2:ntime
    y = X(:,t-1);
    X(:,t)=rk4(dt,y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
end


function dy = couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6)
    dy(1)=-sigma*y(1)+sigma*y(2);
    dy(2)=-y(1)*y(3)+(1+C1*y(4))*k*y(1)-y(2);
    dy(3)=y(1)*y(2)-b*y(3);
    dy(4)=(C2*y(2)+C3*y(5)+C4*y(4)*y(5)-Od*y(4)+Sm+Ss*cos(2*pi*t/Spd))/Om;
    dy(5)=(C5*y(4)+C6*y(4)*y(5)-Od*y(5))/Sigma;

end

function y_next = rk4(dt,y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6)
    dy = couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k1 = dt*dy';
    dy = couple(y+k1/2,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k2 = dt*dy';
    dy = couple(y+k2/2,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k3 = dt*dy';
    dy = couple(y+k3,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k4 = dt*dy';
    y_next = y+k1/6+k2/3+k3/3+k4/6;
end

And I try to use fft to calculate the power spectrum

xn = X(4,:); % X(1,:)X(2,:)X(3,:) for x_1,x_2,x_3, X(4,:) for omega, X(5,:) for eta.
Fs= 1/dt;
nfft = 2 ^ nextpow2(ntime);
cxn = xcorr(xn,'unbiased');  

CXk = fft(cxn);
psd2 = abs(CXk);

index = 0:round(nfft/2-1);
k = index*Fs/nfft;
psd2 = psd2/max(psd2);
psd2 = 10*log10(psd2(index+1));
plot(k,psd2);
xlim([0 100])
grid on

but the results seem different with the reference(https://www.nonlin-processes-geophys.net/24/681/2017/npg-24-681-2017.pdf). On the other hand I don't quite understand how to acquire the confidence upper limit as the reference Figure3 shown. Hope that anyone can give me some advice. Thanks in advance.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.