How can I solve an ode where one of the variables is a matrix?
3 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Are Mjaavatten
2021-7-8
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 个评论
Are Mjaavatten
2021-7-9
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!