Solving system of differential equations
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.
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.