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

조회 수: 2 (최근 30일)
Ikechi Ndamati 2021년 7월 7일
댓글: Ikechi Ndamati 2021년 7월 9일
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 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이전 댓글 1개 표시이전 댓글 1개 숨기기
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.
Ikechi Ndamati 2021년 7월 9일
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.

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

### Community Treasure Hunt

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

Start Hunting!

Translated by