3
$\begingroup$

I am encountering two issues while running the following code in Mathematica:

Slow execution: I have been running the code for approximately 30 minutes, but it is still executing. I am using Windows 64-bit with Mathematica version 13.2.1.0. Is there any way to optimize this code to make it run faster? Are there any specific techniques or alternative functions that I can utilize to expedite the execution time?

Closed form solutions: Additionally, I would like to inquire if it is possible to obtain closed form solutions for the functions m[t], n[t], p[t], and q[t] within the given set of equations.

equations = {
   m'[t] == 1/a (r3 + n[t] - q[t] - (r1 + r2) m[t]),
   n'[t] == -(1/b) m[t],
   p'[t] == 1/c (m[t] - r4),
   q[t] == p[t] + (m[t] - r4) r5,
   m[0] == m0,(*Initial condition for m(t)*)
   n[0] == n0,(*Initial condition for n(t)*)
   p[0] == p0   (*Initial condition for p(t)*)};

sol = DSolve[equations, {m, n, p, q}, t]
Improve this question
$\endgroup$

2 Answers 2

Reset to default
4
+100
$\begingroup$

The system as presented to DSolve consists of autonomous linear DEs plus an algebraic equation that on the surface depend on each other. Consequently, DSolve seems unable to take advantage of the simplicity of the autonomous linear system and applies a DAE solver instead.

Here's an approach: Eliminate the algebraic equation for q[t], solve the linear DE system, and construct the solution for q from the solution to the linear system and the formula for q:

equations = {
   m'[t] == 1/a (r3 + n[t] - (p[t] + (m[t] - r4) r5 ) - (r1 + r2) m[t]),
   n'[t] == -(1/b) m[t],
   p'[t] == 1/c (m[t] - r4),
   
   m[0] == m0,  (*Initial condition for m(t)*)
   n[0] == n0,  (*Initial condition for n(t)*)
   p[0] == p0   (*Initial condition for p(t)*)};

sol = Append[#, q -> Function @@ {t, p[t] + (m[t] - r4) r5 /. #}] & /@
   DSolve[equations, {m, n, p}, t] // AbsoluteTiming
(* {69.7276, <1.7 MB output>} *)

Quite a bit faster than the 30+ min. reported in the OP.

Update 1: A faster DSolve[]

We can reduce the number of parameters and simplify the formulas. This makes a marked improved, much more than I would have expected.

{bmat, amat} = CoefficientArrays[
   First@SolveValues[equations[[;; 3]], D[Through[{m, n, p}[t]], t]],
   Through[{m, n, p}[t]]
   ];

idx = 0;
reducedAmat = 
  Replace[Normal@amat,
   x_ /; ! PossibleZeroQ[x] :> $a[++idx], {2}];
backsub = Simplify@
  DeleteCases[
   Thread[Flatten@reducedAmat -> Flatten@Normal@amat], 
   HoldPattern[0 -> 0]];
(* simplification (by inspection): *)
reducedAmat = reducedAmat /. $a[3] -> -$a[2];

(*ClearSystemCache[]*)(* clear caches for valid timing *)
DSolve[
 {D[{m[t], n[t], p[t]}, t] == 
      reducedAmat . {m[t], n[t], p[t]} + bmat, 
  equations[[4 ;;]]},
 {m[t], n[t], p[t]}, t] /. 
   backsub // AbsoluteTiming

(*  {2.46386, < 0.9 MB output >}  *)

Update 2: Comment on MatrixExp[] solutions

Aside from lacking working code, the answer from @RolandF has a perhaps subtle mathematical error.

The general matrix exponential solution to $\dot {\bf y}(t) = {\bf A} \mathbin{.} {\bf y}(t) + {\bf b(t)}$ has the form $${\bf y}(t)=\Psi(t) \mathbin{.}\left({\bf y}(0)+\int_0^t \Psi(-t) \mathbin{.} {\bf b}(t) \; dt\right) \,,\quad \Psi(t) = \exp({\bf A}\,t)$$ When $\bf b$ is constant and $\bf A$ is invertible, this is easily simplified to a solution resembling @RolandF's in outline; however, $\bf A$ is not invertible in the OP's linear system. And since one of the basis functions of the solution space to the homogeneous system is constant (eigenvalue $0$) and $\bf b$ is constant, the situation is equivalent to having a repeated eigenvalue. There is no way to avoid the integral (or a calculation equivalent to it, since the integral would be simple were it not for the numerous parameters making superficially complicated).

{bmat, amat} = CoefficientArrays[
   First@SolveValues[equations[[;; 3]], D[Through[{m, n, p}[t]], t]],
   Through[{m, n, p}[t]]
   ];

psi = MatrixExp[amat t] // Simplify; // AbsoluteTiming
(*  {2.79531, Null}  *)

msol = psi . ({m0, n0, p0} + # - (# /. t -> 0) &@
      Integrate[(psi /. t -> -t) . bmat, t]); //
    AbsoluteTiming
(*  {4.77207, Null}  *)
    
msol // ByteCount
(*  158312  *)

(* put in solution form *)
matexpsol = 
  Thread[{m, n, p} -> (Function[t, #] & /@
   Transpose@ msol)];
matexpsol = 
   Append[matexpsol, 
      q -> Function @@ {t, p[t] + (m[t] - r4) r5 /. matexpsol}];

matexpsol // ByteCount
(*  262976  *)

If we use reduceAmat as in Update 1, we reduce the times above to a little less than half, which is about as fast DSolve in Update 1.

If you worry about the use of the fundamental theorem of calculus to compute msol, one can, at some cost of time, compute the definite integral:

msol = psi . ({m0, n0, p0} + 
      Integrate[(psi /. t -> -t) . bmat, {t, 0, t}]); //
    AbsoluteTiming
(*  {25.6736, Null}  *)

The result is equivalent (applying Simplify to the difference of the two yields 0).

Update 3: Alternative conversion to linear system

Another way to convert the DAE system to a linear ODE system is to differentiate the equation for q[t]. DSolve computes the solution in about a minute, about as fast as the first solution, but the solution here is much larger in size.

equations3 = {
   m'[t] == 
    1/a (r3 + n[t] - (p[t] + (m[t] - r4) r5) - (r1 + r2) m[t]), 
   n'[t] == -(1/b) m[t], p'[t] == 1/c (m[t] - r4),
   D[q[t] == p[t] + (m[t] - r4) r5, t],
   m[0] == m0,(*Initial condition for m(t)*)
   n[0] == n0,(*Initial condition for n(t)*)
   p[0] == p0,   (*Initial condition for p(t)*)
   q[t] == p[t] + (m[t] - r4) r5 /. t -> 0};

PrintTemporary@Dynamic@Clock@Infinity;
DSolve[equations3, {m, n, p, q}, t] // ByteCount // AbsoluteTiming

(*  {64.1766, 5816160}  *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ very nice and complete solutions. I'll wait for a bit more time before giving the bounty. $\endgroup$
    – internet
    2 days ago
4
$\begingroup$

The set of equations

equations = 
   {m'[t] == 1/a (r3 + n[t] - q[t] - (r1 + r2) m[t]), 
   n'[t] == -(1/b) m[t], p'[t] == 1/c (m[t] - r4), 
   q[t] == p[t] + (m[t] - r4) r5, 
   m[0] == m0,(*Initial condition for m(t)*)
   n[0] == n0,(*Initial condition for n(t)*)
   p[0] == p0   (*Initial condition for p(t)*)};

boils down to a linear equations of first order for three unkown functions

MatrixForm[ 
equations3 = 
Collect[{m'[t] == 1/a (r3 + n[t] - q[t] - (r1 + r2) m[t]),
n'[t] == -(1/b) m[t],p'[t] == 1/c (m[t] - r4)} /. 
{q -> ( p[#] + (m[#] - r4) r5 &)} //FullSimplify, 
{m[t], n[t], p[t]}]]

In tradtional form

$$\partial_t \left(\begin{array}{c} m \\ n \\ p \\ \end{array}\right) =\left( \begin{array}{ccc} -\frac{\text{r1}+\text{r2}+\text{r5}}{a} & \frac{1}{a} & -\frac{1}{a} \\ -\frac{1}{b} & 0 & 0 \\ \frac{1}{c} & 0 & 0 \\ \end{array} \right).\left( \begin{array}{c} m \\ n \\ p \\ \end{array} \right)+\left( \begin{array}{c} \frac{\text{r3}+\text{r4} \text{r5}}{a} \\ 0 \\ -\frac{\text{r4}}{c} \\ \end{array} \right)$$

Quite generally, any non-homogenous linear equation with constant coefficients has as solution the constant solution with x' ==0 plus the general solution of the homogenous equation, it is the exponential of the constant factor. that can be a matrix if x is a vector

   DSolve[x'[t] == a  x[t] + b, x[t], t]

    {{x[t]->-(b/a)+E^(a t) Subscript[\[ConstantC], 1]}}

The Eigensystem of the matrix yields the three (two) frequencies, that multiply the time variable t in the exponentials. The solution has three free constants C[1],C[2],C[3] to be adapted to fit the start conditions.

The expression are indefinitely long by a blow up of the different forms of the square roots in the Eigenvalue. Try to replace all the expressions multiplying the time variable t by a frequency $\pm \omega$, there is only one, because one eigenvalue is 0.

The starting point is

   dsol = DSolveValue[equations, {m[t], n[t], p[t]}, t];

The structure of the solutions can be inspected by setting two constants to zero.

 dsol1 /. {C[2] :> 0, C[3] :> 0};

The frequencies can be determined by

Cases[dsol1, t*__, \[Infinity]] // Flatten//Simplify // Union

Alternatively the solutions can be determined in the direct way

 x' = A x + b  ;  xspec= Solve[ A x == b,x],   xgen = MatrixExp[A (t-t0]] . x[t0]
Improve this answer
$\endgroup$
1
  • $\begingroup$ +1 and thanks for the nice idea. $\endgroup$
    – internet
    2 days ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.