2 views (last 30 days)

Show older comments

Hello, Please I am trying to solve the ODE in the screenshot as part of a larger code that applies the split-step fourier method to the nonlinear schrodinger's equation. The other parts of my code works well except the ODE shown in the screenshot. I need to calculate Nc from the first equation so that I can ultimately calculate my alpha_f. However, Nc is a function of the field A (which is represented in the code as u). However, Nc at any point in time is a function of the field A at that particular time. For example, when z is fixed (each loop of the split-step algorithm is at a particular z) at the 5th time interval, Nc(5) = rho*|u(5)|^4 - b.*N(5). So in other words, I will like the ode to extract the fifth element of u and use it to slove Nc(5)

The initial value of Nc is 0.

I have tried solving the ODE using both ode45 and the runge-kutta (RK4) method without success (you will notice that I commented the ode45 method that I tried). Can you help me point out what I am not doing properly? If I can solve this using any method, I will be glad.

Ps: This code will give error as it is if you run it. This is because of the ODE portion, however, the split-step fourier algorithm works perfectly as I have tested it in isolation.

clear

%% Basic Parameters

%This portion of the code creates the time, frequency and wavelength

%domain/grid for the simulation

%For uniformity of units, all lengths are in km and time are in ps

NT = 2^11; %Number of Pixels(grid points)

lambda_c = 1550e-12; %center of wavelength[km] 1550nm

c = 299792458e-15; %Speed of Ligtht[km/ps]

omega_c = 2*pi*c/lambda_c; %center of angular freq[THz]

vo = c/lambda_c; %central frequency

%% Silicon waveguide parameters

sigma = 1.45e-27; %free carrier coefficient [km^2]

h = 6.63e-52; %Planck's const [km^2*kg/ps]

a_eff = 0.3e6; %effective area [ps^2]

beta_tpa = 5e-15; %TPA coefficient [km/W]

rho = beta_tpa/(2*h*vo*a_eff^2);

b = 1/3e3; %1/3e3

n2 = 6e-24; %Kerr coefficient [km^2/W]

%% Field/grid parameters

t = -0.1:1/NT:0.1;

n = length(t);

f = NT*(-n/2:n/2 - 1)/n; %define bin/freq

omega = (2*pi).*f; %Angular frequency axis

lambda_axis = 2*pi*c./(omega+omega_c); %Wavelength axis

%% Pulse parameters

Po = 25;%3 % peak power of input [W]

t_fwhm = 30e-3; %2 % duration of input [ps]

to = t_fwhm/(2*sqrt(2*log(2)));

Ao = sqrt(Po); %Amplitude

pulse = input('Enter 1 for sech pulse, 2 for gaussian pulse'); % Pulse

%% Fiber parameters

fiber_length = 0.092; %[km]

fiber_loss = 1;%[dB/km]

gamma = ((2*pi*n2)/(lambda_c*a_eff)) +1i*beta_tpa/a_eff; % nonlinear coefficient [1/W/km]

TR = 0.005; %[ps]

alpha_l = fiber_loss/4.343; % attenuation coefficient

fiber_dispersion = -0.0185;%[ps/nm/km]

fiber_d_slope = 0.02;%[ps*ps/nm/km]

beta2 = -0.18e3; %ps^2/km

beta3 = 4; %ps^3/km

Lnl = 1/(gamma*Po)

Ld = to^2/abs(beta2)

tic

%% Pulse field profile

A = Ao*exp(-(t.^2/(2*to^2))); %Gaussian Pulse

%% Convert Pulse to frequency domain

Af = fftshift(fft(A));%Fourier transform of the pulse envelope in freq domain

Af_inp = abs(Af); % Input pulse spectrum

%% Split-step fourier algorithm

loop = 1; %number of loops

n_steps = 200; %number of steps

delta_z = fiber_length/n_steps; %step size

L = (1i*0.5*beta2*omega.^2)+(1i*0.1667*beta3*omega.^3)-0.5*(alpha_l); %Linear operator

%% SSF loop

%Beginning of SSFM loop

for i = 1:loop*n_steps

%step1

%% 1st symmetric linear portion

%Solve linear part in freq domain

A1 = Af.*exp(L*delta_z/2);

%Solve nonlinear part in time domain

u = ifft(ifftshift(A1));

%Beginning of ODE calculation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% RK4

tt = -0.1:1/NT:0.1;

Z = abs(interp1(t,u,tt));

% Z = Z';

AA = @(tt) Z;

t(1) = -0.1;

y(1) = 0;

f = @(t,y) rho.*(abs(AA(t))).^4 + b.*y;

h = 1/NT;

for i = 1:length(t)

%update x

t(i)= tt(i);

t(i+1) = t(i) + h;

%update y

k1 = f(t(i) ,y(i));

k2 = f(t(i)+0.5*h, y(i)+0.5*k1*h);

k3 = f(t(i)+0.5*h, y(i)+0.5*k2*h);

k4 = f(t(i)+h, y(i)+k3*h);

y(i+1,:) = y(i)+(h/6)*(k1+ 2*k2 +2*k3 +k4);

end

%% ODE45

% tt = -0.1:1/NT:0.1;

% Z = abs(interp1(t,u,tt));

% Z = Z';

% AA = @(tt) Z;

% y = zeros (410,1);

% f = @(t,y) rho.*(abs(AA(t))).^4 + b.*y; % Are

% [tt,y] = ode45(f,tt,0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%End of ODE calculation

alpha_f = y*sigma;

%% Nonlinear portion

% N = 1i*gamma*abs(u).^2;

abs2A = abs(u).^2;

absA1 = fftshift(fft(abs(u).^2));

absA2 = fftshift(fft((abs(u).^2).*u));

second_term = (1i/omega_c).*((ifft(ifftshift(-1i*omega.*fftshift(fft(abs2A)))))+conj(u).*ifft(ifftshift(-1i*omega.*fftshift(fft(u)))));

third_term = TR.*ifft(ifftshift(-1i*omega.*absA1));

N = 1i*gamma*(abs2A + second_term - third_term);

%% 2nd Symmetric linear portion

A2 = u.*exp(alpha_f*delta_z);

A2 = A2.*exp(N*delta_z);

v = fftshift(fft(A2));

Af = v.*exp(L*delta_z/2);

Aft = ifft(ifftshift(Af));

Af = fftshift(fft(Aft));

end

% End of SSF loop

Ab = abs(Af);

%% Output plots

%Plot input pulse

figure

subplot (1,2,1)

plot(t,abs(A));

title('Input pulse');

xlabel('time (t)');

ylabel('Amplitude');

%Plot input spectrum

subplot (1,2,2)

plot(lambda_axis,Af_inp/max(Af_inp));

xlim ([1.2e-9 2e-9]);

title('Input Spectrum');

xlabel('Wavelength');

ylabel('Amplitude');

figure

plot(lambda_axis,Af_inp/max(Af_inp),'--');

xlim ([1.2e-9 2e-9]);

hold on

plot(lambda_axis,Ab/max(Af_inp));

xlim ([1.2e-9 2e-9]);

title('Output Spectrum');

xlabel('Wavelength');

ylabel('Amplitude');

figure

plot(lambda_axis,Ab);

xlim ([1.2e-9 2e-9]);

title('Output Spectrum');

xlabel('Wavelength');

ylabel('Amplitude');

figure

imagesc(lambda_axis,i*(1:loop*n_steps)/fiber_length,abs(Ab).^2);

colormap hot

ylabel('L/L_D')

xlabel('T/T_0')

axis square

toc

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!