5 views (last 30 days)

Show older comments

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)

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

Start Hunting!