2
$\begingroup$

I have the stochastic system which consists of 4 nonlinear equations. White Gaussian noise is used in the third equation only. Nevertheless, the whole system is stochastic.

Some problems arise when I try to compute random trajectories:

General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation.

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[…, _SystemException].

There is plenty of free memory, but for some reasons MATHEMATICA considers my system as an uncomputable one. Is there a way to overcome it? here is my code:

(*parameters*)
ρ = 0.2; νmax = 100; Vth = 25; kV = 20; gKLeakOvergl = 0.5; \
GsynOvergl = 5; Kbath = 8.5; τK = 100;  γ = 10; δK \
= 0.02; τNa = 20;  δNa = 0.03; τm = 0.01; τD = \
2; σOvergl = 25; δxD = 0.01;

(*additional functions*)
    Ipump[K_, Na_] = ρ/((1 + Exp[3.5 - K]) (1 + Exp[(25 - Na)/3]));
    f[x_?NumericQ] = 
      Dot @@ MapAt[Boole, 
        Internal`FromPiecewise[Piecewise[{{0, x < 0}, {x, x >= 0}}]], 1];
    ν[V_] := νmax f[2/(1 + Exp[-2 (V - Vth)/kV]) - 1];
    Vk[K_] = 26.6 Log[K/130];
    uOvergL[K_, V_, xD_] = 
      gKLeakOvergl (Vk[K] - Vk[K0]) + GsynOvergl ν[V] (xD - 0.5); 

(*initial conditions*)
K0 = 3;  Na0 = 10; V0 = Vk[K0]; xD0 = 1;

(*differential equations*)
rhs1[K_, Na_, V_, xD_] := (Kbath - K)/τK - 
  2 γ Ipump[K, Na] + δK ν[V]; 
rhs2[K_, Na_, V_, xD_] := (Na0 - Na)/τNa - 
  3 Ipump[K, Na] + δNa ν[V];
rhs3[K_, Na_, V_, xD_] := (-V + uOvergL[K, V, xD])/τm;
rhs4[K_, Na_, V_, xD_] := (1 - xD)/τD - δxD xD ν[V];

proc = ItoProcess[{\[DifferentialD]K[t] == 
    rhs1[K[t], Na[t], V[t], 
      xD[t]] \[DifferentialD]t, \[DifferentialD]Na[t] == 
    rhs2[K[t], Na[t], V[t], 
      xD[t]] \[DifferentialD]t, \[DifferentialD]V[t] == 
    rhs3[K[t], Na[t], V[t], 
       xD[t]] \[DifferentialD]t + σOvergl/τm \
\[DifferentialD]w[t], \[DifferentialD]xD[t] == 
    rhs4[K[t], Na[t], V[t], xD[t]] \[DifferentialD]t}, {K[t], Na[t], 
   V[t], xD[t]}, {{K, Na, V, xD}, {K0, Na0, V0, xD0}}, 
  t, {w \[Distributed] WienerProcess[]}]

solproc = 
  RandomFunction[proc, {0., 100., 0.1}, Method -> "EulerMaruyama"];
$\endgroup$
1
+50
$\begingroup$

First, I changed K to k, since K is a built-in symbol. Then I ran the process for a shorter time and found that V[t] got huge:

solproc = RandomFunction[proc, {0., 1.2, 0.1}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]

Mathematica graphics

I then tried solving the deterministic part with NDSolve, which worked fine:

sol = NDSolve[{
  k'[t] == rhs1[k[t], Na[t], V[t], xD[t]], 
  Na'[t] == rhs2[k[t], Na[t], V[t], xD[t]], 
  V'[t] == rhs3[k[t], Na[t], V[t], xD[t]], 
  xD'[t] == rhs4[k[t], Na[t], V[t], xD[t]], k[0] == K0, 
  Na[0] == Na0, V[0] == V0, xD[0] == xD0}, {k, Na, V, xD}, {t, 0, 10}];

Plot[Evaluate[{k[t], Na[t], V[t], dX[t]} /. sol], {t, 0, 10}]

Mathematica graphics

Going back to the stochastic model, decreasing the step size seems to fix the problem.

solproc = RandomFunction[proc, {0., 100, 0.01}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]

Mathematica graphics

No apparent improvements with a smaller step size:

solproc = RandomFunction[proc, {0., 100, 0.001}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]

Mathematica graphics

$\endgroup$

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.