Substitutions with ODE23/45 in MATLAB

bflbfl

Hi everybody,

I am trying to solve laser rate equations with MATLAB using ODE23/45. For that I have to write the derivative terms always on the left hand side. This creates a lot of redundant variables in the function. Here is a part of the code:

dy(1) = I/(q*V) - y(1)/Tau_n - y(4)*y(2);
% derivative of the no. of carriers, N(t)

dy(2) = (y(4)-1/Tau_p+log(100)/Tau_D) * y(2) + R_sp;
% derivative of the no. of photons, S(t)

dy(3) = (alpha/2)*3.2e3*(y(1)-N_th) + 20/Tau_D;
% derivative of the optical phase, phy(t)

dy(4) = y(5)*(-Kappa*y(6))+(1-Kappa*y(2))*y(7);
% derivative of g(N,S)

dy(5) = 3.2e3 * y(7);
% derivative of g(N)

dy(2) = y(6);
dy(1) = y(7);

Because of y(6) and y(7), I have to define initial values for them as well, for instance y(6) = y(7) = 1. This is now making problems for y(1) and y(2) since they are changing in each iteration, but because of the constant values of y(6) = y(7) = 1, MATLAB considers y(1) = y(2) = 1. Obviously, I do not want this to happen. I cannot write the dy(2) and dy(1) terms on the right hand side, since it is not allowed by MATLAB. Any suggestions to overcome this problem would be welcome. Thanks :)

bflbfl

Just to mention that Tau, Kapa, Tau_n, Tau_p and alpha are constants.

AlephZero

Can you post the mathematical equations you are trying to solve, not just your failed attempt at writing Matlab code?

I am guessing you have some basic misunderstanding about what you need to do here.

bflbfl

Thank you for your reply. So here are the equations:

dN/dt = I/q - N/Tau_n - g(N,S)*S
% derivative of the no. of carriers, N(t)

dS/dt = [g(N,S) - 1/Tau_p + log(100)/Tau_D] * S + R_sp;
% derivative of the no. of photons, S(t)

dphy/dt = (alpha/2) * 3.2e3 * (N-N_th) + 20/Tau_D;
% derivative of the optical phase, phy(t)

dg(N,S)/dt = dg*[-Kappa*g(N)]+[1-Kappa*S]*dS/dt;
% derivative of g(N,S)

dg(N)/dt = 3.2e3 * dN/dt;
% derivative of g(N)

There are more equations in fact but I have written the basic simplified ones, for the sake of understanding.

According to MATLAB rules, one cannot write derivative terms on the right hand side. So I am using substitutions, as per the information available on mathworks website. This is what I did with y(6) and y(7) in the code for dS/dt and dN/dt in the code respectively. I look forward to your suggestions in this connection !!!

AlephZero

I don't see why you need to introduce y(6) and y(7) at all.

If you wanted, you could substutute for dN/dt in
dg(N)/dt = 3.2e3 * dN/dt
to get
dg(N)/dt = 3.2e3 * [ I/q - N/Tau_n - g(N,S)*S ]
and similarly for dg(N,S)/dt.

That "gets rid of the derivatives on the right hand side" without introducing any more variables.

That's the way I would code it in matlab. Just write

dy(1) = I/(q*V) - y(1)/Tau_n - y(4)*y(2);
dy(2) = (y(4)-1/Tau_p+log(100)/Tau_D) * y(2) + R_sp;
dy(3) = (alpha/2)*3.2e3*(y(1)-N_th) + 20/Tau_D;
dy(4) = y(5)*(-Kappa*y(6))+(1-Kappa*y(2))*dy(1);
dy(5) = 3.2e3 * dy(2);

There is no problem using the value of dy(1) and dy(2) again after you have calculated them.

Or if you feel worried about doing that, you could write

temp1 = I/(q*V) - y(1)/Tau_n - y(4)*y(2);
temp2 = (y(4)-1/Tau_p+log(100)/Tau_D) * y(2) + R_sp;
dy(1) = temp1;
dy(2) = temp2;
dy(3) = (alpha/2)*3.2e3*(y(1)-N_th) + 20/Tau_D;
dy(4) = y(5)*(-Kappa*y(6))+(1-Kappa*y(2))*temp1;
dy(5) = 3.2e3 * temp2;