solving time differential equations in frequency domain

If I have a set of time differential equations for rho_ij and I need to find let's say [d(rho_12)/d(omega)], where omega is frequency. Note that these differential equations have terms given as function of omega. what is the best way to do this task. Is it true to use (fft)

[Mathieu NOE]

hello

do you have a document to share with the DE ?

[Mariam]

Dear Mathieu,

Thank you for your reply.

I need to know is it correct the method I have used to get rho_12 as function of freuency and the derivative of rho_12 with respect to the frequency by using the set of time differential equations. I neet your help to determine the best tool to solve this problem. Can FFT be useful towards achieving my aim.

[Mariam]

%-------------------------------------------------------------------------%
function MTMA1_GV(theta_M,theta_G,R,R_M,Omega_c,I_p,Difference)
tic
% The values for general constants in the system:
el_c=1.602e-19;hbar=(6.63e-34)/(2*pi);gamma_M=1e14;%(Damping in MPN)
epsilon_0=8.85e-12; k_B=1.381e-23; E_F=1.36*el_c;c=3e8;T=300;
lambda=570e-9; omegaR=2*pi*c/lambda;% the value of freq. at resonance
omega_pL=1.36e16;% For Plasmon frequency
V_F=1e8;mu=1e4; N=1e20;
%omega_c=omega; omega_p=omega;
%-------------------------------------------------------------------------%
%                                                                         %
%                         The system parameters                           %
%                                                                         %
%-------------------------------------------------------------------------%
% For the dimensions in the system;
L_x=0.5e-9; d=L_x; L_z=7e-9; L_y=L_z;
V_G=L_x*L_y*L_z; %V_G=pi*L_z^2*L_x; 
V_M=(4/3)*pi*R_M^3;
%------------------------------
% For the angles in the system:
phi_1=pi-0.5*pi-theta_G; phi_2=pi-0.5*pi-theta_M;
theta_Q=phi_1+phi_2;
%------------------------------
% For the QD to be used in the system:
mu_12=1e-28;mu_13=mu_12;
gamma_12=1e9; gamma_13=gamma_12; gamma_32=0.35*gamma_12;
gamma_q=0.5*(gamma_12+gamma_13+gamma_32);
gamma=(el_c*V_F^2)/(mu*E_F);%(Damping in Graphene)
%-------------------------------------------------------------------------%
t=0:400e-15:400e-15;% Time interval
epsilon_b=12.9; epsilon_q=6.5; epsilon_inf=5.7;
E_p=sqrt(2*I_p/(epsilon_0*epsilon_q*c));
Omega_p=mu_12*E_p/hbar;
i=0;
%Difference=0.1e15; % the value of small change away from the resonance
Step=((omegaR+Difference*omegaR)-(omegaR-Difference*omegaR))/1000;
for omega_new=omegaR-Difference*omegaR:Step:omegaR+Difference*omegaR
    i=i+1;
    omega=omega_new;
%-------------------------------------------------------------------------%
epsilon_star=(2*epsilon_b+epsilon_q)/(3*epsilon_b);
epsilon_G= 1 + ...
       (1i/(epsilon_0*d*omega))*(el_c^2/(8*hbar)* ...
       (tanh((hbar*omega + 2*E_F)/(4*k_B*T))  + ...
        tanh((hbar*omega - 2*E_F)/(4*k_B*T))) + ...
       (-1i*el_c^2/(8*pi*hbar))* ...
       (log((hbar*omega+2*E_F)^2/((hbar*omega-2*E_F)^2+(2*k_B*T)^2)))+...
       (1i*el_c^2/(pi*hbar))*(E_F/(hbar*omega + ...
       1i*hbar*gamma)));
epsilon_M=epsilon_inf-((omega_pL)^2/((omega)^2+1i*omega*gamma_M));
%-------------------------------
% For the system polarizability:
zeta_x=1 - pi*L_x/(2*L_z);zeta_z=pi*L_x/(4*L_z); %zeta_y= zeta_z;
alpha_Gx=4*pi*V_G*(epsilon_G-epsilon_b)/(3*epsilon_b + ...
                                         3*zeta_x*(epsilon_G-epsilon_b));
alpha_Gz=4*pi*V_G*(epsilon_G-epsilon_b)/(3*epsilon_b + ...
                                         3*zeta_z*(epsilon_G-epsilon_b));
alpha_M=V_M*((epsilon_M-epsilon_b)/(epsilon_M+2*epsilon_b));
%-------------------------------
% For the distances in the system:
R_GM=L_x+L_z+R;
R_QG=(sin(theta_M)/sin(theta_Q)*R_GM);
R_QM=(sin(theta_G)/sin(theta_Q)*R_GM);
%--------------------------------
pi_x=(1/(4*pi*epsilon_star))*(alpha_Gx*(3*cos(phi_1)-1)/R_QG^3 +...
                               alpha_M* (3*cos(phi_2)-1)/R_QM^3); 
pi_z=(1/(4*pi*epsilon_star))*(alpha_Gz*(3*cos(theta_G)-1)/R_QG^3 +...
                               alpha_M* (3*cos(theta_M)-1)/R_QM^3);
%-------------------------------
phi_x=(-1*alpha_Gx*alpha_M/((R_GM^3)*(4*pi*epsilon_star)^2))*...
      ((3*cos(phi_1)-1)/R_QG^3 + (3*cos(phi_2)-1)/R_QM^3);
phi_z=(2*alpha_Gz*alpha_M/((R_GM^3)*(4*pi*epsilon_star)^2))*...
      ((3*cos(theta_G)-1)/R_QG^3 + 3*cos(theta_M)-1/R_QM^3);
%-------------------------------
Lambda_x=(mu_12^2/(2*pi*(4*pi*epsilon_star)^2*(hbar*epsilon_0*epsilon_b)))*...
(alpha_Gx*(3*cos(phi_1)-1)^2/R_QG^6 + alpha_M*(3*cos(phi_2)-1)^2/R_QM^6);
Lambda_z=(mu_13^2/(2*pi*(4*pi*epsilon_star)^2*(hbar*epsilon_0*epsilon_b)))*...
(alpha_Gz*(3*cos(theta_G)-1)^2/R_QG^6 + alpha_M*(3*cos(theta_M)-1)^2/R_QM^6);
%-------------------------------------------------------------------------%
%                                                                         %
%                             The QDs parameters                          %
%                                                                         %
%-------------------------------------------------------------------------%
% For the detuning and coupling constants:
omega_12=2.172436*(el_c/hbar);
omega_13=2.172432*(el_c/hbar);
%-------------------------------------------------------------------------%
% In this code we don't need delta_p & delta_c
omega_p=omega; %omega_c=omega;
Delta_p=omega_12-omega_p;
Delta_c=0;%omega_13-omega_c;
%-------------------------------------------------------------------------%
%-------------------------------------------------------------------------%
rho0=[0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00]; % The initial values 
%[rho,fval]=fsolve(@Stability,rand(9)./10,options);% rand to generate random 
                                                   % values b/w 0-1
[t,rho]=ode45(@f,t,rho0); % The function to find the sol for diff eqs.
rho12=interp1(t,(rho(:,2)),400e-15); % To get sol for rho12 in freq.
rho11=interp1(t,(rho(:,1)),400e-15);
rho22=interp1(t,(rho(:,5)),400e-15);
PopInv(i)=rho11-rho22;
Chi(i)=(N*mu_12/(epsilon_0*E_p))*rho12;
x_axis(i)=Delta_p; % For x axis being detuning
ResonanceFreq(i)=omega;  % For x axis being frequency 
P(i)=(rho12);
end
D=diff(P);% Differential of P=rho12;
VG=c./(1+0.5.*ResonanceFreq(:,1:length(D)).*real(D));% Group Velocity
nG=(c./VG)-1;% Group Index
%-------------------------------------------------------------------------%
%save DisPlusAbs1.mat
plot(x_axis(:,1:length(D))./1e12,nG,'blu')
%plot(x_axis./1e12,PopInv,'bla');
%plot(x_axis./1e12,imag(Chi),'bla');%hold on
%plot(x_axis./1e12,real(Chi),'blu');hold off
yline(0)
xlabel('\Delta_p THz','fontsize',12)
%ylabel('Absorption Im [\chi]','fontsize',12)
%legend('Absorption','Dispersion')
%plot(x_axis./1e15,VG./c,'-bla');hold on
%plot(x_axis./1e15,nG,'-blu');hold off
%set(gca,'fontsize',12)
%-------------------------------------------------------------------------%
% In case for plotting in both sides of y axis with one x axis:
%colororder({'bla','blu'})
%yyaxis left
%plot(x_axis./1e15,VG./c,'-bla');hold on
%ylabel('c/v_g')
%yyaxis right
%plot(x_axis./1e15,nG,'-blu');hold off
%ylabel('n_g')
%xlabel('The frequency \omega GHz','fontsize',12)
%xlabel('\Delta_p GHz','fontsize',12)
set(gca,'fontsize',12)
%-------------------------------------------------------------------------%
%                                                                         %
%                        Equations of motion for the system               %
%                                                                         %
%-------------------------------------------------------------------------%
% Note that: rho11 --> rho(1), rho12 --> rho(2), rho13 --> rho(3),
%            rho21 --> rho(4), rho22 --> rho(5), rho23 --> rho(6),
%            rho31 --> rho(7), rho32 --> rho(8), rho33 --> rho(9),
%-------------------------------------------------------------------------%
    function MEs = f(t,rho)
MEs=zeros(9,1);
% The eqns of motion:
MEs(1)= -(gamma_12+gamma_13)*rho(1)+1i*Omega_c*(pi_z+ phi_z)*rho(7)+...
                                     1i*(Lambda_z)*rho(3)*rho(7)+...
                                     1i*Omega_p*(pi_x+phi_x)*rho(4)+...
                                     1i*(Lambda_x)*rho(2)*rho(4)-...
                                     1i*Omega_c*conj(pi_z+phi_z)*rho(3)-...
                                     1i*conj(Lambda_z)*rho(7)*rho(3)-...
                                     1i*Omega_p*conj(pi_x+phi_x)*rho(2)-...
                                     1i*conj(Lambda_x)*rho(4)*rho(2);
MEs(2) = -(0.5*(gamma_13)+0.5*(gamma_12))*rho(2)-1i*(Delta_p)*rho(2)+1i*(Lambda_x)*...
                                     (rho(5)- rho(1))*rho(2)+...
                                     1i*(Omega_p)*(pi_x+phi_x)*(rho(5)-rho(1))+...
                                     1i*(Omega_c)*(pi_z+phi_z)*rho(8)+...
                                     1i*(Lambda_z)*rho(3)*rho(8);
MEs(3) = -(0.5*(gamma_13)+0.5*(gamma_12)+0.5*(gamma_32))*rho(3)-...
                                     1i*(Delta_c)*rho(3)+...
                                     1i*(Lambda_z)*(rho(9)- rho(1))*rho(3)+...
                                     1i*(Omega_c)*(pi_z+phi_z)*(rho(9)-rho(1))+...
                                     1i*(Omega_p)*(pi_x+phi_x)*rho(6)+...
                                     1i*(Lambda_x)*rho(2)*rho(6);
MEs(4) = -(0.5*(gamma_13)+0.5*(gamma_12))*rho(4)+1i*(Delta_p)*rho(4)-1i*conj((Lambda_x))*...
                                    (rho(5)- rho(1))*rho(4)-...
                                     1i*(Omega_p)*conj(pi_x+phi_x)*(rho(5)-rho(1))-...
                                     1i*(Omega_c)*conj(pi_z+phi_z)*rho(6)-...
                                     1i*conj((Lambda_z))*rho(7)*rho(6);
MEs(5) = (gamma_12)*rho(1)+(gamma_32)*(rho(9)-rho(5))-1i*(Omega_p)*(pi_x+phi_x)*rho(4)-...
                                     1i*(Lambda_x)*rho(2)*rho(4)+...
                                     1i*(Omega_p)*conj(pi_x+phi_x)*rho(2)+...
                                     1i*conj((Lambda_x))*rho(4)*rho(2);
MEs(6) = -(0.5*(gamma_32))*rho(6)+1i*(Delta_p-Delta_c)*rho(6)-1i*(Omega_c)*(pi_z+phi_z)*rho(4)-...
                                     1i*(Lambda_z)*rho(3)*rho(4)+...
                                     1i*(Omega_p)*conj(pi_x+phi_x)*rho(3)+...
                                     1i*conj((Lambda_x))*rho(4)*rho(3);
MEs(7) = -(0.5*(gamma_13)+0.5*(gamma_12)+0.5*(gamma_32))*rho(7)+...
                                     1i*(Delta_c)*rho(7)-...
                                     1i*conj((Lambda_z))*(rho(9)-rho(1))*rho(7)-...
                                     1i*(Omega_c)*conj(pi_z+phi_z)*(rho(9)-rho(1))-...
                                     1i*(Omega_p)*conj(pi_x+phi_x)*rho(8)-...
                                     1i*conj(Lambda_x)*rho(4)*rho(8);
MEs(8) = -(0.5*(gamma_32))*rho(8)-1i*(Delta_p-Delta_c)*rho(8)+1i*(Omega_c)*conj(pi_z+phi_z)*rho(2)+...
                                     1i*conj((Lambda_z)*rho(7)*rho(2)-...
                                     1i*(Omega_p)*(pi_x+phi_x)*rho(7)-...
                                     1i*(Lambda_x)*rho(2)*rho(7));
MEs(9) = (gamma_13)*rho(1)+(gamma_32)*(rho(5)-rho(9))-1i*(Omega_c)*(pi_z+phi_z)*rho(7)-...
                                     1i*(Lambda_z)*rho(3)*rho(7)+...
                                     1i*((Omega_c)*conj(pi_z+phi_z))*rho(3)+...
                                     1i*conj((Lambda_z))*rho(7)*rho(3);
%MEs(10) = rho(1) + rho(5) + rho(9) -1;
    end
%-------------------------------------------------------------------------%
toc
end