5
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$ – Geoff Oxberry Jun 29 '13 at 20:35
  • $\begingroup$ Hi. deltat is a time step. And in the model the exponential factor for x'3 is e^((t/deltat)^2)) $\endgroup$ – CarlosCr Jun 29 '13 at 22:45
  • 1
    $\begingroup$ Having the time step for a numerical method inside your ODE system makes no sense; that information is for the discrete formulation only. $\endgroup$ – Geoff Oxberry Jun 30 '13 at 19:48
5
$\begingroup$

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$.

Improve this answer
$\endgroup$
12
  • $\begingroup$ 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. $\endgroup$ – CarlosCr Jun 29 '13 at 22:50
  • $\begingroup$ @user4671 You have to have initial condition, initial condition means what the system is like when you start your timer. $\endgroup$ – Shuhao Cao Jun 29 '13 at 23:21
  • $\begingroup$ Depending on the argument to the exponential, ode45 may be inappropriate due to stiffness, and replaced with ode15s. $\endgroup$ – Geoff Oxberry Jun 30 '13 at 20:14
  • $\begingroup$ The matlab code that i wrote for this problem is the next: i´ve included some constants for simulation purpouses: $\endgroup$ – CarlosCr Jun 30 '13 at 21:03
  • 1
    $\begingroup$ en.m.wikipedia.org/wiki/Nondimensionalization $\endgroup$ – Geoff Oxberry Jul 1 '13 at 2:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.