Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Substitutions with ODE23/45 in MATLAB

  1. Oct 5, 2011 #1
    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 :)
     
  2. jcsd
  3. Oct 5, 2011 #2
    Just to mention that Tau, Kapa, Tau_n, Tau_p and alpha are constants.
     
  4. Oct 5, 2011 #3

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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.
     
  5. Oct 6, 2011 #4
    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 !!!
     
  6. Oct 6, 2011 #5

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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;
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Substitutions with ODE23/45 in MATLAB
Loading...