MATLAB Answers

How to solve an ODE of a function of a function

2 views (last 30 days)
Ikechi Ndamati
Ikechi Ndamati on 14 Jul 2021
Commented: Ikechi Ndamati on 18 Jul 2021
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.
%% 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)
%% 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
%% 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);
%% 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 of SSF loop
Ab = abs(Af);
%% Output plots
%Plot input pulse
subplot (1,2,1)
title('Input pulse');
xlabel('time (t)');
%Plot input spectrum
subplot (1,2,2)
xlim ([1.2e-9 2e-9]);
title('Input Spectrum');
xlim ([1.2e-9 2e-9]);
hold on
xlim ([1.2e-9 2e-9]);
title('Output Spectrum');
xlim ([1.2e-9 2e-9]);
title('Output Spectrum');
colormap hot
axis square
Ikechi Ndamati
Ikechi Ndamati on 18 Jul 2021
Thanks @Torsten I was able to eliminate the errors and get it to run by changing AA = interp1(t_axis,u,t); to AA = interp1(u,t);
Now I have an output (which I think is incorrect. I'd just have to look closely) just by implementing the above. However, when I implemented the suggestion, it did not bring any error. However, there was no error but the output was blank.
Additionally, please how does the ode45 output work? Why is it that the output of my y is a 137x410 matrix instead of 1x410 that I expect it to be?

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!