6

I need to solve the following system of differential equations. enter image description here

When I have the solutions for $n_f$ and $v$, I need to find and plot $J=-e_\cdot n_{f} \cdot v$.

I wrote a code in matlab with all ODEs like this:

      function systemSolve
      clc


      tr=50e-12;        % Recombination lifetime
      n0=1e11;         % Density of free carriers (Recordar que es 1e17 cm-3)
      tc=2e-12;         % Trapping time
      ts=30e-15;        % Carrier scattering time
      m=0.067*9.11e-31; % effective mass GaAs
      ev=8.854e-12;     % permitivity  
      n=900;            % factor geometry
      q=1.6e-19         % electron charge

      de=50e-15         %  delta t


      timeRange=[0 0.1e-12]; 
      initialConditionVector=[0;1e-15;1e-15;1e-15];
      [t,x]=ode45(@xprime,timeRange,initialConditionVector);
      figure(1),plot(t,x(:,1))

      J=q*x(:,3).*x(:,1);
      figure(2),plot(t,J)

      function f=xprime(t,x)
      f=[x(4); ...
        -(x(2)/tr)+(x(3)*q*x(1)); ...
        -(x(3)/tc)+(n0*(exp((t/de)^2))); ...
        -((1/ts)*x(4))-(((x(3)*q^2)/m*ev)*x(1)/n)+(q*x(2)/(m*n*ev*tr))];
      end
end

I suppose that:

$x_1=v$

$x_2=P_{sc}$

$x_3=n_f$

$x_4=x'_1$

I expect to find a current pulse like figure 1 but I get a exponential solution.

Figure 1

What is wrong in this code and can you suggest me how i can solve this problem?

CarlosCr
  • 61
  • 3
  • Welcome to SciComp! Could you elaborate a bit on the physics of your model? Is "deltat" a problem parameter, or does that refer to a time step? – Geoff Oxberry Jun 29 '13 at 20:35
  • Hi. deltat is a time step. And in the model the exponential factor for x'3 is e^((t/deltat)^2)) – CarlosCr Jun 29 '13 at 22:45
  • 1
    Having the time step for a numerical method inside your ODE system makes no sense; that information is for the discrete formulation only. – Geoff Oxberry Jun 30 '13 at 19:48

1 Answers1

5

You can use MATLAB's ode45, please refer to the first example here: Solve nonstiff differential equations; low order method.

Here I use $y$ instead of $x$. If you don't know how to write a $\mathtt{func.m}$ script and save it in your path, you can use function handle as well, in your case:

% output of function handle in vector form
rigid = @(t,y)[x(4); -x(2)+x(3);  ??? ; ???];
options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
[t,y] = ode45(rigid,[0 T],[y10 y20 y30 y40],options);

??? is left for you to fill up, you can tweak the options if you read the document of the link above, y10 through y40 are the initial values. Lastly J = y(:,1).*y(:,3) then plot $J$ against $t$.

Shuhao Cao
  • 2,552
  • 17
  • 30
  • Hi. Thanks a lot for your assistance. I don´t have initial values in my problem. it is possible avoid this initial conditions?. or I have to deduce them. – CarlosCr Jun 29 '13 at 22:50
  • @user4671 You have to have initial condition, initial condition means what the system is like when you start your timer. – Shuhao Cao Jun 29 '13 at 23:21
  • Depending on the argument to the exponential, ode45 may be inappropriate due to stiffness, and replaced with ode15s. – Geoff Oxberry Jun 30 '13 at 20:14
  • The matlab code that i wrote for this problem is the next: i´ve included some constants for simulation purpouses: – CarlosCr Jun 30 '13 at 21:03
  • tr=50e-12;%Recombination lifetime n0=1%e11;%Density carriers tc=2e-12%Trapping time ts=30e-15;%scattering time m=0.0679.11e-31;%mass LT-GaAs ev=8.854e-12;%vaccum n=900;%geometrical factor q=1.6e-19%electron de=50e-15%deltat timeRange=[0 0.1e-12]; initialConditionVector=[0;1e-15;1e-15;1e-15]; [t,x]=ode45(@xprime,timeRange,initialConditionVector); figure(1),plot(t,x(:,1)) J=qx(:,3).x(:,1); figure(2),plot(t,J) function f=xprime(t,x) f=[x(4);-(x(2)/tr)+(x(3)qx(1));-(x(3)/tc)+(n0(exp((t/de)^2)));-((1/ts)x(4))-(((x(3)q^2)/mev)x(1)/n)+(qx(2)/(mnevtr))]; end end – CarlosCr Jun 30 '13 at 21:10
  • With this code I did not found the solution that I hope. What is wrong? – CarlosCr Jun 30 '13 at 21:11
  • @user4671: If possible, you should nondimensionalize your equations so that the rescaled state variables and time variable are roughly 1 over the new time interval of interest. – Geoff Oxberry Jul 01 '13 at 00:08
  • @GeoffOxberry I don't understand your answer. can you help me explain it? – CarlosCr Jul 01 '13 at 02:11
  • @ShuhaoCao I already do it. can you help me? – CarlosCr Jul 01 '13 at 02:11
  • 1
    http://en.m.wikipedia.org/wiki/Nondimensionalization – Geoff Oxberry Jul 01 '13 at 02:25
  • The line n0=1%e11; % Density of free carriers (Rqcordar quq qs 1q17 cm-3) is very peculiar. Is n0 supposed to be 1 or 1e11? – Bill Greene Jul 01 '13 at 16:27
  • @BillGreene Yes it is 1e11 – CarlosCr Jul 01 '13 at 16:44