# Solving system of differential equations

In solving the following system using Mathematica, I get

DSolve::bvfail: For some branches of the general solution, unable to solve the conditions. >>

The equations are

$${dx\over dt}=\lambda -\beta x v-d x$$

$${dy\over dt}=\beta x v-a y$$

$${dv\over dt}=-uv$$
where $$\lambda, \beta, d,a,u$$ are constant. The Mathematica code is

DSolve[{x'[t] == lambda - d*x[t] - beta*x[t]*v[t],    y'[t] == beta*x[t]*v[t] - a*y[t], v'[t] == -u*v[t], x[0] == xstar,    y[0] == ystar, v[0] == vstar}, {x[t], y[t], v[t]}, t]


According to the book I am reading, the output I should get is

Assuming that $$x$$ is constant on the relevant time scale,
$$y(t)=y(0){ue^{-at}-ae^{-ut}\over u-a}$$
$$v(t)=v(0)e^{-ut}$$

Why don't I get a solution and what is the meaning of "Assuming that $$x$$ is constant on the relevant time scale". Does that mean x is at steady state?

• DSolve can't solve many nonlinear systems, because no closed form solutions exist. The statement about assuming $x$ is constant is an extra assumption the authors make to get a solution. More background info would be helpful in interpreting that statement: what's the model of and what is the book? Sep 3, 2017 at 14:39
• If x is actually constant, dx/dt is zero and this is a linear system. But there is information missing -- how can the solution not involve beta or lambda? What constant value is assumed for x? Sep 3, 2017 at 14:45

Update Added going beyond the text perturbatively with comparison against numerics.

# Recovering Text-book

So the problem you're running into is that Mathematica's just not able to solve the differential equations exactly given the constraints you've offered.

Let's first see if we can indeed meet your book's approximation, which does hold x is in a steady state; it's derivative is zero. x[t]=x[0]=xstar.

xEqns = {x'[t] == lambda - d*x[t] - beta*x[t]*v[t],  x[0] == xstar}

vEqns = {v'[t] == -u*v[t],  v[0] == vstar}

yEqns = {y'[t] == beta*x[t]*v[t] - a*y[t],  y[0] == ystar}

constXConds= {x'[t] -> 0, x[t] -> xstar ,  x[0] -> xstar}

constXSol =
DSolve[Flatten[{yEqns, vEqns}] /. constXConds, {y[t], v[t]}, t] //
FullSimplify // Flatten


yields:

{v[t] -> E^(-t u) vstar,
y[t] -> (E^(-t (a + u)) (beta E^(a t) vstar xstar -
E^(t u) (beta vstar xstar + (-a + u) ystar)))/(a - u)}


which looks complicated relative to your book's solution (typos corrected).

Let's see what conditions we need to match to the book's simpler (and one less parameter $\beta$). I don't know why certain parameters (like $\beta$) aren't present in your books solution, so let's see if these two can match at all:

Your book's solution with obvious typo corrected:

bookSoln = {v[t] == vstar Exp[-u t],
y[t] == ystar ( u Exp[-a t] - a Exp[- u t]) / (u - a)}


Matching occurs when the following is zero:

letsSee = bookSoln /. constXSol // FullSimplify


{True, (E^(-t (2 + u)) (E^(2 t) - E^(t u)) (beta vstar xstar - 2 ystar))/(-2 + u) == 0}

We applied our solution constXSol (which was in the form of a rule mapping $y(t)\to$ our solution) to the RHS of the book soln. The resulting equations have to be zero to match. We can just use Solve here to get the constraint on beta (if we weren't able to already read it off by eye). i.e. Solve[letsSee, beta]. The equation here is really not difficult. But let me show you a tool that should be in your toolbox when you've got complicated sometimes transcendental functions in variables:

If two expressions can be equal they will be equal order by order in their series expansion around their independent variable. I.e. if we series expand both our solution and the book solution around $t$, then order by order they must equal. As we only have one free parameter that shows up in our solution compared to the books, if these guys can be equal they'll be equal at the first order in t. So we solve the series expansion of letsSee instead of the full equation.

someSol = Flatten@Solve[Normal@Series[letsSee, {t, 0, 1}], beta]


beta -> (a ystar)/(vstar xstar)}

This is sufficient to be correct to all orders:

Indeed: letsSee /. someSol // ExpandAll $\mapsto$ 0.

Now what's happened here? For some reason your book chose the initial velocity of y (y') to vanish and set Beta accordingly. You can verify by plugging in:

Flatten[yEqns /. someSol /. t -> 0 /. {y[0] -> ystar, v[0] -> vstar, x[0] -> xstar} // FullSimplify] $\mapsto$ {y'[0] == 0, True}

Note: y of course y still evolves even with steady state (x[t]->star) conditions (c.f time dependence of the book solution or constXSol /. someSol // FullSimplify).

# Going Beyond Text

Going beyond the text to first order.

So let's see how to use Mathematica to go a little beyond your textbook now.

First lets get some numbers so we can compare against a numerical integration.

someParams = {xstar -> 1, ystar -> 1, vstar -> 1, u -> 1/101,
a -> 1/99, d -> 1/100, xprime[0] -> 1};


I want to be able to specify the initial velocity of x to be xprime, so I'll solve for lambda in terms of this:

initXCond =
xEqns[[1]] /. t -> 0 /. x'[0] -> xprime[0] /.
x[0] -> xstar /. v[0] -> vstar // Solve[#, lambda] & //
Flatten


{lambda->d xstar+beta vstar xstar+xprime[0]}

Since I personally don't care about lambda, let's just set that forever.

initXCond /. Rule :> Set;


Recall, we needed to set beta to speak to your text solution, so lets do that now

Echo[someSol] /. Rule :> Set;


{beta->(a ystar)/(vstar xstar)}

OK, We're good to solve numerically:

numSoln =
NDSolve[{xEqns, yEqns, vEqns} /. someParams, {x, y, v}, {t, 0,
100}];


First let's see how good an approximation our assume x is constant solution was.

Plot[{{y[t], v[t]} /.
numSoln, {y[t], v[t]} /. x[t] -> xstar /.
Rule @@@ bookSoln /. someParams} // Flatten // Evaluate, {t, 0,
4}, PlotLegends->{"numeric y","book y"}]


For these values doesn't look great. Not so unexpected, I personally chose x'[0] nontrivial given other constants. Good let's now see what going beyond the book gets us.

The fact that we can consider even a constant x for a while means we can assume a polynomial form in t for x[t].

xApprox[n_] :=
Sum[fudge[i] t^i, {i, 0, n}] /. fudge[0] -> xstar /.
fudge[1] -> xprime[0]
xRulex[n_] := x[t] -> xApprox[n];


Lets see that the obvious choice of fudge parameters satisfy the xEqns to 0th order at t=0

xEqns /. Map[D[#, t] &, xRulex[12]] /. xRulex[12] /. t -> 0 /.
v[0] -> vstar /. x[0] -> xstar


{True,True}

Now let's see what we need to satisfy to go one higher order.

First, notice the $v$ equations are separable, so let's solve them:

vOnly = DSolve[vEqns, {v[t]}, t] // Flatten


{v[t]->E^(-t u) vstar}

First order xEqns are given:

Series[xEqns[[1]] /. Map[D[#, t] &, xRulex[2]] /. xRulex[2] /.
vOnly, {t, 0, 1}]


xprime[0]+2 fudge2 t+O[t]^2==xprime[0]+(a u ystar-d xprime[0]-(a ystar xprime[0])/xstar) t+O[t]^2

Let's solve this in order t, which will inform parameters up to fudge2 because of the derivative.

x1eqn = Normal@
Series[xEqns[[1]] /. Map[D[#, t] &, xRulex[2]] /. xRulex[2] /.
vOnly, {t, 0, 1}] /. Equal[a_, b_] :> a - b //
FullSimplify


t (2 fudge2+d xprime[0]+(a ystar (-u xstar+xprime[0]))/xstar)

x1soln = Solve[x1eqn[[2]] == 0, fudge[2]] // Flatten //
FullSimplify


{fudge2->(a u xstar ystar-(d xstar+a ystar) xprime[0])/(2 xstar)}

Now lets go for y assuming 1st order in t for x

y1Soln = DSolve[yEqns /. xRulex[1] /. vOnly, y[t], t] //
FullSimplify // Flatten


{y[t]->1/((a-u)^2 xstar) E^(-a t) ystar (-a u xstar+u^2 xstar+a xprime[0]-a E^(t (a-u)) (xprime[0]-(a-u) (xstar+t xprime[0])))}

Now lets go for y assuming 2nd order in t for x. (We don't need to see what it does in x, it'll just mean an appropriately defined fudge3.)

y2Soln = DSolve[yEqns /. xRulex[2] /. vOnly /. x1soln, y[t], t] //
FullSimplify // Flatten


{y[t]->1/(2 (a-u)^3 xstar^2) E^(-t (a+u)) ystar (-2 E^(t u) u^3 xstar^2+a^4 E^(a t) t^2 ystar (u xstar-xprime[0])+a^3 E^(a t) (2 xstar (xstar-t u (1+t u) ystar)+t ((2-d t) xstar+2 (ystar+t u ystar)) xprime[0])+a xstar (2 E^(t u) (2 u^2 xstar+d xprime[0]-u xprime[0])+E^(a t) (2 u^2 xstar+2 u (1+t u) xprime[0]-d (2+t u (2+t u)) xprime[0]))+a^2 (-2 E^(t u) (xstar+ystar) (u xstar-xprime[0])+E^(a t) (u xstar (-4 xstar+(2+t u (2+t u)) ystar)+(2 (-1+d t+t (-2+d t) u) xstar-(2+t u (2+t u)) ystar) xprime[0])))}

The moment of truth, let's see how we did.

Plot[{(y[t] /.
numSoln), (y[t] /. y1Soln /. xprime[0] -> 0 /. vOnly /.
someParams), (y[t] /. y1Soln /. vOnly /. someParams), (
y[t] /. y2Soln /. vOnly /. someParams)} // Flatten //
Evaluate, {t, 0, 25},
PlotLegends -> {"y - numeric",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, $$0$$]]\) (book)",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, $$1$$]]\)",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, $$2$$]]\)"}]


Looks like we're doing great, but of course we're interested in precision science. Let's see how far we go to get stay within percent of the actual solution:

Plot[Map[(# - (y[t] /. numSoln))/(y[t] /. numSoln) &, {(y[t] /.
numSoln), (y[t] /. y1Soln /. xprime[0] -> 0 /. vOnly /.
someParams), (y[t] /. y1Soln /. vOnly /. someParams), (
y[t] /. y2Soln /. vOnly /. someParams)}] // Flatten //
Evaluate, {t, 0, 25}, PlotRange -> {-.02, .02},
PlotLabel -> "Percentage Deviation", Frame -> True,
PlotLegends -> {"Numeric",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, \
$$0$$]]\)/\!$$\*SubscriptBox[\(y$$, $$num$$]\)",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, \
$$1$$]]\)/\!$$\*SubscriptBox[\(y$$, $$num$$]\)",
"\!$$\*SubscriptBox[\(y$$, SubscriptBox[$$x$$, \
$$2$$]]\)/\!$$\*SubscriptBox[\(y$$, $$num$$]\)"}]


So compare against $|.01|$ to see where our perturbative solution remains within 1% of the numeric integration.

• Thank you for the answer. but I don't understand what is done with the lines letsSee = bookSoln /. constXSol // FullSimplify; and someSol = Flatten@Solve[Normal@Series[letsSee, {t, 0, 1}], beta]. What is {t, 0, 1} in this line and I don't understand what is been done here and why. Would you please be able to clarify this Sep 4, 2017 at 2:12
• As I understand you found $\beta$ which would give the simplified version as in the book.What I did was I set x to be constant at say xstar and tried the following code.DSolve[{y'[t] == beta*xstar*v[t] - a*y[t], v'[t] == -u*v[t], y[0] == ystar, v[0] == vstar}, {y[t], v[t]}, t] The output I get is y[t] -> (E^(-a t - t u) (beta E^(a t) vstar xstar - beta E^(t u) vstar xstar + a E^(t u) ystar - E^(t u) u ystar))/(a - u). If I substitute beta with what you found as {beta->(a ystar)/(vstar xstar)} I get the result in the book. Sep 4, 2017 at 5:49
• By looking at the equations I beleive the {beta->(a ystar)/(vstar xstar)} is obtained under the steady state solution for y. But if x is held constant why isn't beta found by setting $\lambda -\beta x v-d x=0$ Sep 4, 2017 at 5:52
• Hi Sam, you're absolutely right. I should have commented on this, all the steady state equation on x does is relate Beta to Gamma, but your book chooses its parameters (for some presumably physical reason or for fun) to have initial $y'$ vanish for steady state x. This $y$ still eventually evolves because $v$ evolves. Later today I'll update my answer with additional pedagogical notes for the mma games I played you might not be familiar with. Please don't hesitate to request help with specific bits that are confusing. Sep 4, 2017 at 6:09
• @JohnJosephM.Carrasco I am using the system keyword in Dsolve. I define a matrix A = {{-F*sd, 0, 0}, {F*sd, -lt - F*st, F*sh}, {0, lt, -F*sh}}, put this in a system form system = X'[t] == A.X[t];, with X[t_] = {d[t], tr[t], h[t]} and finally solve it sol = DSolve[system, {d, tr, h}, t], however I don't know how to define the initial conditions. Could you help me with that. In principle, if I would be able to define the IC, is it possible to plot the solution? Dec 4, 2019 at 10:52

You can directly solve this system with DSolve, if you split it into two steps, since v-equation can be solved separately.

eqs = {x'[t] == lambda - d*x[t] - beta*x[t]*v[t],
y'[t] == beta*x[t]*v[t] - a*y[t], v'[t] == -u*v[t], x[0] == xstar,
y[0] == ystar, v[0] == vstar};

vsol = v /. First@DSolve[{v'[t] == -u*v[t], v[0] == vstar}, v, t]

(*   Function[{t}, E^(-t u) vstar]   *)

eqs2 = DeleteCases[eqs /. v -> vsol, True]

{xsol[a_, beta_, d_, lambda_, u_, vstar_, xstar_, ystar_],
ysol[a_, beta_, d_, lambda_, u_, vstar_, xstar_, ystar_]} =
{x, y} /. First@DSolve[eqs2, {x, y}, t]


You get solutions for x and y with unevaluated integrals, which are automatically evaluated when inserting numbers for parameters. In the case analytic integral could not be evaluated, you can switch to NIntegrate at certain t.

xsol[1, 1, 1, 1, 1, -1, 1, 1][t]

(*   -E^(1 - E^-t - t) (-1 + (E - E^(1 + 1/E) -
ExpIntegralEi[1] + ExpIntegralEi[1/E])/
E - (-E^(1 + 1/E) + E^(E^-t + t) + ExpIntegralEi[1/E] -
ExpIntegralEi[E^-t])/E)   *)

Plot[Evaluate[xsol[1, 1, 1, 1, 1, -1, 1, 1][t]], {t, 0, 10},
PlotRange -> All]


(xsol[1, 1, 1, 1, 1, -1, 1, 1] /. Integrate -> NIntegrate)[2]

(*   1.37564   *)


Now you can play around to find the solution of the text book.