How can I solve an ode where one of the variables is a matrix?

Please I cant figure out how to solve an ode where one of the variables is a matrix. The ode is of the form:

dy(t)/dt = rho*|A(t)|^4 + tau*y(t); where rho and tau are constants and A(t) is a 1xn matrix.

I have tried both the Runge-Kutta and ode45 but I think I am not representing the equations rightly. Please can anyone advise me on how to rectify it? My code is shown below:

NT = 2^10; %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);

tau = 1/3e3;

%% Field/grid parameters

tspan = -10:1/NT:10;

n = length(tspan);

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

to = 2;

A = sqrt(3)*exp(-(tspan.^2/2*to.^2)); %Gaussian Pulse

% Declaring the ODE

h = 1/NT;

y(1) = 0;

f=@(x,y) rho.* abs(A(x)).^4 - tau*y;

%%RK4

for i = 1:ceil(xfinal/h)

%update x

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

%update y

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

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

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

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

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

end

plot (x,y);

Answer 1

[Are Mjaavatten]

A(t) should be a function, not a vector. Below I define A as an anonymous function. I also use an anonymous function for your f function.

NT = 2^10;                                  %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);
tau = 1/3e3;
to = 2;
A = @(t) sqrt(3)*exp(-(t.^2/2*to^2)); %Gaussian Pulse

f = @(t,y) rho*abs(A(t))^4 + tau*y;
[t,y] = ode45(f,[-10,10],0);
figure;plot(t,y)

[Ikechi Ndamati]

Thank you @Are Mjaavatten. Your answer worked perfectly fine. But I noticed that this will be difficult to apply it in the real code that I am working on, so let me share the full code:

So instead of using A as you suggested, I want to use u, which is actually A after some operations have been carried out on it (Kindly see code).

In the code below, I changed the A function to be A = @(t) u(t) because I actually want to replace that function with u. Sorry that I did not share the code initially. I thought it will be too cumbersome.

Ps: the ODE portion is within the array of '%' signs.

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^10;                                  %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 = -10:1/NT:10;
    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]
to = 28.4e-3; %2             % duration of input [ps]
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 = 16; % nonlinear coefficient [1/W/km]
gamma = ((2*pi*n2)/(lambda_c*a_eff)) +1i*beta_tpa/a_eff;
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; %-(lambda_c^2*fiber_dispersion)/(2*pi*c);%23.6e-3; %ps^2/km
beta3 = 4; %0.02*lambda_c^4/(2*pi*c)^2; %ps^3/km
Lnl = 1/(gamma*Po)
Ld = to^2/abs(beta2)
% alpha_f = sigma*N_c;
tic
%% Pulse field profile
%This enables the selection of either sech or gaussian pulse
if pulse == 1
   A = Ao*sech(t/to); % Sech pulse, input field [W^(1/2)] 
else
%  C = -5;
%  m = 1;
%  A = Ao*exp(-(1+1i*C*0.5)*(t/to).^(2*m)); %Gaussian Pulse
   A = Ao*exp(-(t.^2/2*to.^2)); %Gaussian Pulse
end
%% 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+alpha_f); %Linear operator

%% loop
for i = 1:loop*n_steps
  %step1

%% 1st symmetric linear portion
L = (1i*0.5*beta2*omega.^2)+(1i*0.1667*beta3*omega.^3)-0.5*(alpha_l);
%
    %Solve linear part in freq domain
    A1 = Af.*exp(L*delta_z/2);
    %Solve nonlinear part in time domain
    u = ifft(ifftshift(A1));
    v = u;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ODE
% A = @(t) sqrt(3)*exp(-(t.^2/2*to^2)); %Gaussian Pulse
A = @(t) u(t)

f = @(t,y) rho*abs(B(t))^4 + b*y;
[t,y] = ode45(f,[-10,10],0);
alpha_f = y*sigma;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Nonlinear portion
        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 = u.*exp(N*delta_z);
      v = fftshift(fft(A2));
      Af = v.*exp(L*delta_z/2); 
      Aft = ifft(ifftshift(Af));
      Af = fftshift(fft(Aft));
      
end
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.48e-9 1.62e-9]);
title('Input Spectrum');
xlabel('Wavelength');
ylabel('Amplitude');

[Are Mjaavatten]

NOTE: This is a revised version of my earlier l comment

The current code is too complex and difficult for me to analyse. There are so many Fourier transforms and inverse transforms that I lose track. Are you sure you do not get lost yourself?

I will just mention a few observations.

Vector for pulse == 1 gives a very narrow pulse, while for pulse == 2 the pulse is much wider that the t interval of -10 to 10.

u(t) in line 77 is a complex vector and not a function, so it cannot be used in ode45. You might use

A = @(tt) = abs(interp1(t,u,tt));
f = @(t,y) rho*abs(A(t))^4 + b*y;   % Are
[tt,y] = ode45(f,[-10,10],0);

Note that, since t is already defined, I use another variable (tt) for the result of ode45. Otherwise it may be confusing.

A general advice is to split your code into functions that do a well-defined operation, and test that each function does what you intend.

[Ikechi Ndamati]

Haha! @Are Mjaavatten I have been on this code for months, so I have gotten used to it and hence, I don't get lost anymore. This is still a first draft of the code and I will clean it up when I get it to run perfectly.

Thank you so so much! I will do that. I will try the interp1.