How can I solve an ode where one of the variables is a matrix?
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Ikechi Ndamati
el 7 de Jul. de 2021
Comentada: Ikechi Ndamati
el 9 de Jul. de 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);
0 comentarios
Respuesta aceptada
Are Mjaavatten
el 8 de Jul. de 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)
3 comentarios
Are Mjaavatten
el 9 de Jul. de 2021
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.
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!