How to solve an ODE of a function of a function

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

[Ruchika Mishra]

Hi,

I found this resource which shows how to calculate ODE using Runge-Kutta 4th order method:

Runge-kutta

please check if this helps

[Ikechi Ndamati]

@Ruchika Mishra Thanks for pointing this out. I know about this code already. My problem is not just how to solve ODE using the RK4, but how to solve an ODE where one of the variables is itself a function of time. I wish to solve the equation in the screenshot.

[Ikechi Ndamati]

So the part of my code that I need help in is:

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

[Ikechi Ndamati]

I need to find the ode of f = @(t,y) rho.*(abs(AA(t))).^4 + b.*y; where AA(t) is supposed to be the interpolation of u. In other words, I really want to solve f = @(t,y) rho.*(abs(u(t))).^4 + b.*y; as in the screenshot. But since u itself is a function of time, it is difficult to achieve.

[Ikechi Ndamati]

Or let me explain it this way @Ruchika Mishra:

I want to solve

dy/dt = rho*u^4 -b*y

where u= u(t)

But now, u(t) is already in the form of a matrix (there is no feasible equation to represent it in because it is a field in time that has already passed through different stages) of 1xN matrix where the elements in N represent the different field powers at different time interval. The initial y = 0. So in every loop iteration, to calculate y, I have to use the corresponding u value at that point in time.

So in each iteration, in order to calculate y from dy/dt, I have to use the previous value of y and the value of u at that point in time.

Recall the equation:

dy/dt = rho*u^4 -b*y

[Torsten]

Use Matlab's ODE45 and look at the example "ODE with time-dependent terms" in the documentation of the solver.

[Ikechi Ndamati]

Thanks @Torsten I have actually done that. However, I was faced with several errors that I could not debug because most of the errors were within the ode45 function. Also, my functions v and u became NaN and the order is now 41x410 instead of 1x410. To be honest, I have no clue on where the problem comes from. I have tried a lot of trial and error unsuccessfully. Please I will appreciate a third eye as I am still a novice that took on a project that is becoming more than I could chew. These are the error messages - I have googled every single one of them, but no success yet.

Error using interp1>reshapeAndSortXandV (line 447)

LENGTH(X) and SIZE(V,1) must be the same.

Error in interp1 (line 128)

[X,V,orig_size_v] = reshapeAndSortXandV(X,V);

Error in SiWv_bugger>myode (line 206)

AA = interp1(t_axis,u,t);

Error in SiWv_bugger>@(t,y)myode(t_axis,y,u,t,rho,b) (line 126)

f = @(t,y) myode(t_axis,y,u,t,rho,b);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in SiWv_bugger (line 127)

[t,y] = ode45(f, tspan, y0, opts);

 t = -0.1:1/NT:0.1;
 tspan = [1 410];
 y0 = 0;%zeros(1,n);
 opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
 f = @(t,y) myode(t_axis,y,u,t,rho,b);
 [t,y] = ode45(f, tspan, y0, opts);
 
function dydt = myode(t_axis,y,u,t,rho,b)
AA = interp1(t_axis,u,t); 
dydt = rho.*(abs(AA)).^4 + b.*y;
%dydt = dydt(:);
end

[Torsten]

Why do you set tspan = [1 410] if you can only provide u in the range [-0.1 : 0.1] ?

[Ikechi Ndamati]

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?