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

3 views (last 30 days)
Ikechi Ndamati on 7 Jul 2021
Commented: Ikechi Ndamati on 9 Jul 2021
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);

Are Mjaavatten on 8 Jul 2021
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 on 9 Jul 2021
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.