40 Finding a “solution”

As you saw in the previous chapter, each differential equation in a dynamical system relates the derivative of a function to the function itself. For instance, in \[\partial_t x = x\,(1-x)\] left-hand side of the equation involves the function $x(t)$. The equation dictates that whatever $x(t)$ might be, it has to be such that $\partial_t x(t)$ is exactly equal to the function $x(t)\,\left(\strut 1- x(t)\right)$ Solving a differential equation is the phrase used for finding such self-consistent functions. This chapter is about techniques for finding solutions.

In high-school mathematics, many algebraic problems were of the form, “Solve for $x$.” In such problems, the answer you got was all you were looking for: the complete answer to the problem.

Having been exposed to this style in high school, you may have come to think of “finding a solution” as the only task to perform in working on a problem.

In working with dynamical systems, “finding a solution” is one of the several tasks, among others, that you might need to perform to serve the purpose of your modeling effort. But not every modeling purpose requires “finding a solution,” that is, finding the function $x(t)$. For instance, the goal of working with a dynamical system is often to sort out other properties of the system, such as the existence, location, and stability of fixed points. (we will introduce fixed points in Chapter 41.)

Given what you learned in studying Block 3 of this book, you likely will be tempted to approach the task of “finding a solution” by applying symbolic anti-differentiation techniques to the problem. After all, each differential equation in a dynamical system involves a function $\partial_t x(t)$. To find $x(t)$ from $\partial_t x(t)$ seems like a matter of applying the “fundamental theorem of calculus,” namely

\[\int \partial_t x(t) dt = x(t)\ .\] Following this logic, we would translate the equation to \[x(t) = \int x(t)\, (1-x(t))dt\ .\] But the problems in Block 3 were all of the form $\frac{d\color{magenta}{x}}{d\color{blue}{t}} = g(\color{blue}{t})$, whereas the problems we work with in this Block are generally of the entirely different form $\frac{d\color{magenta}{x}}{d\color{blue}{t}} = g(\color{magenta}{x})$. Thus, we will usually need special techniques suited to the format of dynamical systems.

Admittedly, in mathematics it is common to refer to integrating a differential equation, but this should be broadly understood as accumulating the increments $\partial_t x(t)$ starting at some initial condition $x(t_0)$, even if that accumulation is not carried out by symbolic anti-differentiation.

In this chapter we will introduce three different techniques to accumulating a solution to a differential equation or a pair of such equations. First, we will look again at the graphical method of “following the flow” in a plot of the flow field. This technique is mainly of use for developing intuition about the dynamics.

Second, we will develop a simple Euler method for accumulating a solution. Third, we will explore how to take a guess about the solution and, when the guess is good enough, refine that into an actual solution. This is called the method of ansätze.

Third, and briefly, we will look at some of the situations where symbolic anti-differentiation can be used. This includes a very brief introduction to substitution methods.

A differential equation like \[\partial_t x(t) = x(t)\,(1-x(t))\] is very busy typographically. One reason is the repeated $(t)$ which play no role other than to state explicitly that $x$ is a function of $t$. A convenient shorthand simply replaces $x(t)$ to make the “is a function of $t$” implicit. The shorthand form of the equation appears as \[\partial_t x = x (1-x)\ .\] it is the reader’s responsibility, knowing that she is working with a differential equation, to remember that $x$ is a function of $t$.

In many texts and research papers, differential equations are written using Leibniz’s notation, like this:

\[\frac{dx}{dt} = x (1-x)\ .\]

This is exactly equivalent to our notation $\partial_t x$.

An even more concise notation, originated by Isaac Newton, is to replace the $\partial_t$ with a simple dot over the function being differentiated as in $\Large\dot{x}$. With this notation the equation looks like \[\dot{x} = x (1-x)\ ,\] about as simple as it gets.

This dot notation is even more expressive when working with second-order differential equations involving second derivative. In the dot notation, $\partial_{tt} x$ is written $\ddot{x}$. Dot notation is also in wide use. Again, it is exactly equivalent to our $\partial_t x$.

40.1 The flow field

With a pair of differential equations, as with the pendulum or the rabbit-fox model, each equation gives one component of the change in state. To draw the flow at single point in state space, evaluate the dynamical functions at that point. Each dynamical function contributes, as its output, one of the components of the state velocity vector. If the parameters in the model have been assigned numerical values, the result of evaluating the right-hand sides will be two numbers.

A case in point is the rabbit-fox system. The axes in the rabbit-fox state space are denominated in units of rabbit density $r$ and fox density $f$. The differential equations are \[\begin{eqnarray} \partial_t r & = & 0.66 r - 1.33 r f\\ \partial_t f & = & -f + rf\\ \end{eqnarray}\]

To to find the state velocity at, say, $r=2, f=1/4$, plug those values into the right-hand side:

$\partial_t r = 1.33 - 0.66 = 0.66\ \ \ $rabbit density per month

$\partial_t f = -0.25 + 0.5 = 0.25\ \ \ $fox density per month.

Once you know the numerical vector value of the state velocity, you need to convert it to a form suitable for plotting in the state space. The conversion is needed because the state space is denominated in rabbit density and fox density, not in rabbit density per month or fox density per month. The conversion is accomplished by multiplying the state velocity vector by a small $dt$, say, 0.1 months.

The conversion produces a vector whose components are denominated in the same way as the state space and thus can be plotted meaningfully in the state space.

To illustrate, let’s draw a flow vector for the state space coordinate $(r=2, f = 1/4)$. Above, we already calculated the components of the state velocity vector;Given the value $\partial_t f = 0.25$ and $\partial_t f = 0.66$. For the sake of illustration, we will set $dt = 0.1$ month. Consequently, the vector to be plotted will be $(0.25, 0.66) dt = (0.025, 0.066)$$ with units of rabbit density and fox density respectively. the right. This flow arrow is drawn in Figure 40.1.

To draw the entire flow field, repeat this process at many other points in the state space as in Figure 40.2.

(a) The flow field depicted by drawing the state velocity vector (multiplied by $dt = 0.1$ to turn it into a length) for many points in the state space.

Some people prefer a visualization of short segments of actual trajectories, as in the right panel in Figure 40.2, rather than the state velocity vector. This is a matter of personal preference.

With the flow field depicted in sufficient detail, you can now trace out trajectory.

To trace out a trajectory, select a initial condition for the system. Then follow the flow, taking only a small step in state space. The next step should be in the direction of the flow arrow at the end of the previous step.

The trajectory you draw will be only a sketch, but it can be effective for developing intuition. Figure 40.3 shows a semi-automated version of the go-with-the-flow method. The computer has been used to draw the arrows. When you click in the plot, the computer also undertakes calculation of the trajectory.

Figure 40.3: The flow field for the rabbit/fox dynamics. Click at an initial state to generate the trajectory from that state. You may need to pinch in or out to see the flow arrows clearly.

Regrettably, from such a sketch of the trajectory, you cannot easily construct $r(t)$ and $f(t)$ for time-series plots. Also, you don’t get a sense of how slow or fast the flow is going. Click at different initial conditions in the flow and you will see different trajectories, each of which is a closed loop, the sort of cycles seen in the dice-free Chutes and Ladders game. But the shape of the trajectory does not tell you whether it takes a long time or a short time to complete a loop.

The next section will show you how the computer constructed the trajectory and how we can get information on the speed of the flow.

40.2 Euler method

Recall from Block 2 the limit definition of the derivative: \[\partial_t x(t) = \lim_{dt \rightarrow 0} \frac{x(t + dt) - x(t)}{dt}\ .\] we will use this definition to develop a very general way to solve differential equations: the Euler method.

The differential equations specify the values of $\partial_t x(t)$ in terms of the dynamical function. In Block 2, we paid attention to whether the limit exists. But here, we know it must because the dynamical functions themselves don’t involve limits. In working with the differential equation it suffices to pick some small, finite $dt$. How small? Pick $dt$ to be small enough that the result wouldn’t change in any substantial way if we used an even smaller time increment, say $dt/10$.

Our starting point for solving each differential equation is to re-write it as a finite difference. To illustrate, we will solve the equation $\partial_t x = x (1 - x)$, which is often called the logistic equation.

Applying the finite difference definition, we get

\[\underbrace{\frac{f(t + dt)- f(t)}{dt}}_{\text{finite-difference approx.}} = \underbrace{x (1-x)}_{\text{dynamical function}}\ .\]

Multiplying both sides of the above by $dt$ and re-arranging terms produces

\[\underbrace{f(t + dt)}_{\text{future state}} = \underbrace{f(t)}_{\text{current state}} +\ \ \ \underbrace{x (1-x) dt}_{\text{step}}\] We call this last equation the Euler formula.

To use this, we start at the initial condition, say $x(t=0) = 0.2$. This initial condition gives us the first row of a tabular representation of the function $x(t)$:

time	state
0	0.2

Next, pick a value for $dt$ that we will use for all the following steps, each of which will add a new row to the table. For the example, we will set $dt = 0.1$. When we have constructed the whole table we can go back and check whether that was small enough.

To fill in the next row, we apply the Euler formula. Sine $dt = 0.1$, the next time step will be $0.1$. Plug in the current state—which is 0.2 right now—to calculate the future state. The step will be $0.2 (1-0.2)\, dt = \color{brown}{0.016}$. Add this step to the current state to get the future state. The table now looks like this:

time	state
0.0	$0.2$
0.1	$0.2 + \color{brown}{0.016} = \color{blue}{0.216}$

The next step will bring us to time $0.2$. Use the Euler formula, pluggin in the value of the present state, $\color{blue}{0.216}$, to find the step. Here that will be $0.216 (1-0.216)\, dt = \color{magenta}{0.0169.}$. Now the table looks like

time	state
0.0	$0.2$
0.1	$0.2 + \color{brown}{0.016} = \color{blue}{0.216}$
0.2	$\color{blue}{0.216} + \color{magenta}{0.0169} = 0.2329$

Add as many rows to the table as you like; the process will be the same.

You will recognize this as an iterative process, as discussed in Chapter 32.

$\ $

As is so often the case, it is wise to think about carrying out processes in terms of fundamental tasks accomplished by calculus operations—evaluate, differentiate, anti-differentiate, solve, find argmax, iterate. The obvious choice for integrating differential equations is “anti-differentiate,” but as described previously, the techniques we covered in Block 3 are not sufficient for the task. Instead, we use iteration to solve differential equations.

This example uses the software you have already seen, Iterate(), to carry out the task. In practice, however, you will use a special form of Iterate() called integrateODE() that makes use of interpolation techniques to give a more precise answer.

To implement the iteration to solve $\partial_t x = x (1-x)$, we need to create a function that takes the current state as input and produces the next state as output. Our one-step function can be this: ::: {.cell layout-align=“center” fig.showtext=‘false’}

next_step <- function(t, x, dt=0.1) {
  t <- t + dt
  x <- x + x*(1-x)*dt
  
  c(t=t, x=x) # return value
}

Notice that we wrote next_step() with an input slot for $dt$. This will not be part of the state being iterated, just a parameter that allows us easily to explore different values for $dt$.

Use Iterate() to carry out the iteration of next_step(). Note that we use the fargs argument to Iterate() to pass our selected value for dt to the function next_step(). We will run the iteration for 100 steps. With $dt=0.1$, those 100 steps will 10 units of time.

Soln <- Iterate(next_step, x0=c(t=0, x=0.2), n=100,
                fargs=list(dt=0.1))

n	t	x
0	0.0	0.2000000
1	0.1	0.2160000
2	0.2	0.2329344
... 101 rows in total ...
99	9.9	0.9998595
100	10.0	0.9998736

We can now plot the time series $x$ vs $t$:

Figure 40.4: The time series by the Euler method with $dt=0.01$.

:::

In the previous example using Iterate() to solve a differential equation, the output of the iteration was a data frame containing values for the solution at discrete times: 0, 0.1, 0.2, and so on. A data table is a perfectly good way to represent a function, but it is handier to have a function in a form that operations like slice_plot() and D() can be applied to. Another way to look at things is that, mathematically, the solution to a differential equation should be a continuous-time function. Fortunately, we have at hand the interpolation techniques covered in Chapter 33 to carry out the construction of a continuous-time function from a tabular representation. The R/mosaic function integrateODE() connects together the iteration and interpolation to provide a solution that is in the form of continuous-time function(s).

Use the R/mosaic function integrateODE() to solve differential equations numerically. It is a specialized function that handles sets of first-order differential equations, but any high-order differential equation can be separated into a set of first-order equations.

To illustrate, this command will solve the differential equation $\partial_t x = x (1-x)$ that we took on in the previous example with Iterate(). ::: {.cell layout-align=“center” fig.showtext=‘false’}

Soln2 <- integrateODE(dx ~ x*(1-x), x = 0.2,
                      bounds(t=0:10), dt=0.01)
## Solution containing functions x(t).

The first argument is a tilde expression, but in a form that is different from from that used in functions such as D() or contour_plot(), etc. To the left of the tilde is a single name composed of the state variable—x here—prefixed by a d. The d is just a reminder that we are describing not x itself, but $\partial_t\ \mathtt{x}$. On the right of the tilde is the function from the differential equation, in this case, $x(1-x)$.

The next argument is the initial condition. We are starting the integration at $x=0.2$. The bounds() sets the time interval for the integration and dt sets the time step..

The output of integrateODE() is an R structure of a type called a “list” that is new to us. The list contains the function(s) created by integrateODE() which you refer to by name (x) using a special form of R punctuation $ suited to lists.. In other words, Soln2$x will be a function, which you can plot like any other function, for instance:

slice_plot(Soln2$x(t) ~ t, bounds(t=0:10))

This use of the R symbol $ is new to us. We won’t emphasize it here. Instead, we’ll use the traj_plot() graphics function (introduced in Chapter 48) which already knows how to access the functions created by integrateODE().

traj_plot(x(t) ~ t, Soln2, bounds(t=0:10))

An important feature of integrateODE() is its ability to handle sets of first-order differential equations. For instance, the rabbit/fox system \[\partial_t r = 0.66\, r - 1.33\, r f\\ \partial_t f = -f + rf\] will be integrated by this command: ::: {.cell layout-align=“center” fig.showtext=‘false’}

Eco_soln <- integrateODE(
  dr ~ 0.66*r - 1.33*r*f, 
  df ~     -f +      r*f,
  r = 2, f = 0.5, #initial conditions
  bounds(t=0:5), dt=0.1)
## Solution containing functions r(t), f(t).

::: You can plot the time series using slice_plot()

pA <- traj_plot(r(t) ~ t, Eco_soln, bounds(t=0:5)) %>% gf_labs(title="Rabbits")
pB <- traj_plot(f(t) ~ t, Eco_soln, bounds(t=0:5)) %>% gf_labs(title="Foxes")
gridExtra::grid.arrange(pA, pB, ncol=2)

Figure 40.5: Times series of the rabbit and fox densities.

To plot the trajectory, simply change the tilde expression used in traj_plot(). which creates a time series plot, traj_plot() shows the trajectory.

traj_plot(f(t) ~ r(t), Eco_soln, nt=10)

:::

40.3 Symbolic solutions

Occasionally it is possible to integrate a differential equation using symbolic techniques. This is particularly true for differential equations that are linear. The example we will handle here is the first-order linear differential equation \[\partial_t x = a\, x\ .\] An advantage of symbolic solutions is that parameters can be handled symbolically.

A method we will use repeatedly in this block is called the “method of ansätze.” An ansatz (singular of the German “ansätze”) is, in this context, a guess for the solution. Since differential equations have been a central part of science for more than 200 years, you can imagine that a large library of equations and their solutions has been assembled. For the equations that are most frequently used and that can be solved symbolically, the solutions are already known. Thus, the “guess” for the solution can be a very well informed guess.

Let’s see how this works for $\partial_t x = a\, x$. From experience, the ansatz will be an exponential function of time, which we can write $x(t) \equiv A e^{\omega t}$. We don’t yet know what is the value of $\omega$ or $A$, so we plug the ansatz into both the left and right sides of the differential equation to work it out.

Plugging in the ansatz, translates the differential equation to a new form:

\[\underbrace{A \omega e^{\omega t}}_{\partial_t x(t)}\ =\ \underbrace{a A e^{\omega t}}_{a x(t)}\ .\] Cancelling out the terms that appear on both sides of the equation gives \[\omega = a\ \ \ \text{which implies}\ \ \ x(t) = A e^{a t}\ .\]

The ansatz substitution didn’t give any result at all for $A$. That is to say, unlike $\omega$, the $A$ is not determined by the differential equation itself. This means that $A$ must be related to the initial condition. Setting $t=0$ gives $x(0) = A$, so in this simple differential equation, $A$ is the initial condition.

A slightly more complex differential equation is \[\partial_t x = a\, (x - b)\ .\] This also has an exponential solution. It is easiest to see this by defining a new variable $y \equiv x - b$. By the rules of differentiation, $\partial_t y = \partial_t x$, so the differential equation can be re-written in the form $\partial_t y = a y$. We already know the solution to this is $y(t) = y_0 e^{a t}$. Translating by to $x$ we get

\[x(t) - b = (x_0 -b) e^{at}\ \ \ \implies x(t) = (x_0 -b)\,e^{at} + b)\ .\]

For nonlinear dynamical function, there is no perfectly general way to find symbolic solutions. But for some dynamical functions, it can be done. we will demonstrate by integrating $\partial_t x = x (1-x)$. The method is made more plausible by using the Leibnizian notation for derivatives, with which the differential equation has this form: \[\frac{dx}{dt} = x(1-x)\ .\] The Leibnizian notation can be interpreted as the ratio of two differentials: $dx$ and $dt$ in this case.

The idea of separating the differential equation is to algebraically move all the $x$ terms to the left side of the equation and all the $t$ terms to the right and then to integrate each side of the equation.

\[dx = x(1-x) dt \ \ \ \implies \ \ \ \frac{1}{x(x-1)}dx = dt\ \ \ \implies\ \ \ \int\frac{1}{x(x-1)}dx = \int dt .\]

The integral on the right side, $\int dt$, should be easily recognizable, giving $\int dt = t + F$, where $F$ is the “constant of integration.”

The integral on the left side may not be as familiar, but the person solving this problem for the second time will remember that \[\frac{1}{x(1-x)} = \frac{1}{x} + \frac{1}{1-x}\] as you can confirm by putting the right side over a common denominator. Each of $1/x$ and $1/(1-x)$ have integrals that are logs: $\int dx/x = \ln(x) + D$ and $\int dx/(1-x) = - \ln(1-x) + E$. Putting the equation back together again, produces

\[\ln(x) + D - \ln(1-x) + E = t + F\ .\]

At this point, move all the constants of integration over to the right side and consolate them into a single constant of integration $C$. At the same time, collect together the two logarithmic terms, giving: \[\ln\left(\frac{x}{1-x}\right) = t + C\ .\] Exponentiate both sides to get:

\[\frac{x}{1-x} = \underbrace{e^C}_{A} e^t\ .\] Since $e^C$ is just a constant, we will write it more simply as $A$.

Now we have \[x = Ae^t - x A e^t \ \ \implies\ \ \ x (1 + Ae^t) = Ae^t\] which gives our solution \[x = \frac{Ae^t}{1 + Ae^t}\ .\] To find the initial condition symbolically, plug in $t=0$, giving $x_0 = A/(1+A)$ or, equivalently $A = x_0/(1-x_0)$. Our previous examples used $x_0 = 0.2$, for which $A = 0.2/0.8 = 0.25$. Graphing this solution gives us the familiar sigmoid:

Symb_soln = makeFun(A*exp(t)/(1 + A*exp(t)) ~ t)
slice_plot(Symb_soln(t, A=0.2) ~ t, bounds(t=-5:10))

Not all differential equations can be separated in this way, and even for those that can, the integrals may not be tractable. So this route to a solution is not a general-purpose one, unlike the Euler method. Still, the Euler method gives only an approximate solution, so with Euler we need to take care that the approximation is close enough for the purpose at hand. In this case, we have both an Euler solution (with $dt=0.1$) and a symbolic solution. Figure 40.6 shows the difference between the two solutions, which ideally should be zero. To show more of the time domain of the solution, we will reset the initial condition to $x_0 = 0.01$. This corresponds to $A = 1/99$.

Figure 40.6: The difference between the Euler and the symbolic solution to $\partial_t x = x (1-x)$ as a fraction of the symbolic solution. At the worst, the Euler solution is off by 1.5 parts in one-million.

40.4 Exercises

Exercise 40.01

Some of these are first-order differential equations, others are not. For each, say whether it is indeed a first-order differential equation and, if so, what is the name of the state variable.

$\partial_t u = f(t)$
$\partial_x y = g(x)$
$\partial_t z = h(z)$
$\partial_x y = f_2(y)$
$\partial_t t = 1$
$\partial_t u = g_2(u) + g_3(u)$

Exercise 40.03

For the following, you can use whatever you like for function names, e.g. $f_1(), f_x(), g(), h(), h_3()$ and so on. Indicate the names of function inputs by writing them explicitly inside the parentheses in the normal way, e.g. $g_2(u, v, w)$.

Write down the framework for a system of differential equations with state variables $u, v, w$.
Write down the framework for a system of differential equations with state $z, u$.
Write down the framework for a system of differential equations with state variables $v, x, z$.

Exercise 40.05

Figure 40.5 shows time series for the rabbit and fox population density starting at the initial condition $r=2, f=0.5$ for the first 5 time units thereafter. The fox graph looks like a hump function, the rabbit graph show a little uptick near $t=5$.

Using the R/mosaic commands given in the text to make Figure 40.5, integrate the equations from $t=0$ to $t=15$ and plot the time series for both rabbits and foxes.

A. Using the time-series plots, estimate the period of the cyclic oscillations. - What is the period of the fox population cycle? - How large in amplitude (peak to trough) is the fox population cycle? - How do the cycle period and amplitude for the rabbits compare to those for the foxes?

B. Change the initial condition from $r=2, f=0.5$ to $r=1, f=0.5$ and plot the time series. - What is the period of the fox population cycle? - How large in amplitude (peak-to-trough) is the fox population cycle?

Exercise 40.06

We will explore some common pitfalls for using the integrateODE() and traj_plot() functions.

Part A What is the purpose of the integrateODE() function?

To translate a dynamical system into a set of numerical solutions, one for each state variable.
To draw a flow field of the dynamical system.
To anti-differentiate a differential equation, giving a symbolic solution.
To plot the trajectory of a dynamical system.

Part B What is the proper format to use in specifying the dynamical system to integrateODE()?

A set of tilde expressions, each one an individual argument.
A set of named arguments, with the formula on the right-hand side.
A set of equations separated by semi-colons.
A set of tilde expressions separated by semi-colons.

Part C What should the left-hand side of the tilde expression look like for a state variable z?

dz z zdot dt_z

Part D For the dynamical system \[\partial_t x = y\\\partial_t y = -x ,\] what should the integrateODE() command look like? (Note: we are using ... as a placeholder for other inputs that will be needed by integrateODE().)

integrateODE(dx ~ y, dy ~ -x, ...)
integrateODE(dx = y, dy = -x, ...)
integrateODE(dx ~ y, dy = -x, ...)

Part E Aside from the tilde-expressions encoding the dynamical system, what other information is absolutely essential to specify when constructing a numerical solution?

The initial condition for each and every state variable.
The trajectory.
The time series.
The initial condition for at least one of the state variables.

Part F How should you specify how long you want integrateODE() to carry forward the solution in time, say for 10 time units?

bounds(t=0:10) duration=10 t=10 for=10

Exercise 40.07

This system of differential equations, called the Lorenz equations is famous as the first widely recognized example of a type of motion called chaos.

What are the state variables and what are the parameters?
What is the dimension of the state space?
For initial conditions $x_0 = -5.276, y_0 = -7.152, z_0=19.452$ and $\rho=28$, $\sigma=10$, and $\beta=8/3$, use integrateODE() to integrate the equations over $0 \leq t \leq 100$ with a stepsize of $dt=0.01$. (This means you should set use the arguments bounds(t=0:50), dt=0.01. Call your trajectory T1.

## Solution containing functions x(t), y(t), z(t).
## Solution containing functions x(t), y(t), z(t).

Based on your results in (3), what are the values of $x$, $y$, and $z$ at time $t=50$.
Make a time series plot of $x(t)$. Note that $x(t)$ jumps between a oscillation with negative $x$ and the same kind of oscillation with positive $x$.
Plot out the $y$ component of the trajectory versus the $x$ component of the trajectory. To get a smooth plot, you will have to use the npts=10000 argument to traj_plot(). The trajectory will appear to cross itself, but this is because you are only plotting two of the state variables.
Create another trajectory T2 in the same way you made T1. But change each of the initial conditions in the third decimal point. Then plot out the $x(t)$ time series for T1 (with color="blue") and the $x(t)$ time series for T2 (with color=“red”`) on the same axes. The two time series will track one another very closely, but then come out of sync with one another.
The characteristic feature of chaos is that a small change in initial conditions can lead to a very large change in the time series. How long does it take for the two time series in (7) to become utterly out of sync with one another?
The function traj_plot_3D() provides a simple way to look at the trajectory in 3 dimensions. Such a plot shows that the trajectory does not in fact cross itself. Use this command:

traj_plot_3D(x, y, z, T1, npts=5000)

Exercise 40.09

The differential equation $\partial_t x = a x$ is so commonly encountered that you should simply memorize the solution: \[x(t) = A e^{a t}\] which you can recognize as exponential growth (or decay, if $a < 0$) from an initial condition $x_0 = A$.

Exponential growth is considered fast, but there are far faster forms of growth. To illustrate, consider the differential equation \[\partial_t x = a x^2\ .\]

This can be interpreted in terms of a model of the size of the flame as one lights a match. Think of the flame as a ball of hot gas of radius $x$; the gas include oxygen (O_2_), nitrogen, and carbon dioxide, as well as a vapor of the combustible material such as the potassium chlorate of the match head.

Within the ball of flame, O_2_ reacts with the cumbustible material to produce the products of combustion and heat. Needless to say, this reaction eliminates the O_2_ in the ball. But O_2_ can diffuse into the ball from outside. The O_2_ infusion rate available in this way is proportional to the surface area of the ball, that is, $a x^2$. Thus the differential equation models the growth of the flame ball.

The match-flame equation is one that can be separated into parts: all the $x$ components on one side, the $t$ on the other. That is:

\[\underbrace{\frac{dx}{dt} = a x^2}_{\text{original DE}}\ \ \ \implies \underbrace{\frac{dx}{x^2} = a dt}_{\text{separated DE}}\]

Integrating both sides of the separated equation will produce $a t + C$ on the right side.

A. Integrate the left side of the separated equation and use that to find a relationship between $x$ and $a t + C$.

B. The constant of integration, $C$, will reflect the initial condition. Plug in $t=0$ and calculate from the relationship in (A) what is $C$ in terms of $x_0$.

C. Replace $C$ in your relationship with the expression in terms of $x_0$ you found in (B). Confirm that this is \[x(t) = \frac{1}{a t - 1/x_0}\ .\]

D. $x(t)$ has a vertical asymptote. Find it.

E. Use integrateODE() to integrate the original differential equation. You will have to pick some numerical value for $a$ and $x_0$. Take care to make $dt$ small enough. You will know that $dt$ is small enough when you get the approximately same solution using $dt/10$.

F. Describe in everyday words what the solution says and how big the ball of flame becomes.

The model $\partial_t x = a x^2$, like all models, includes some features of the real-world system and excludes others. In this case, once the ball reaches a critical diameter, there is no longer enough combustion product to continue the reaction at the rate depicted in the model. If you watch a match being lit, you will see both the explosion and the eventual exhaustion of the combustion material.

Exercise 40.11

This activity makes use of the following app:

Click on the picture of the app and it will open in a new browser tab. Arrange that new tab side-by-side with the one where you are reading this.

To solve a differential equation with the Euler method, you need two things:

The differential equation itself. Several choices are available in the selector on the left side of the app. On the right side of the equation is the dynamics(x) function.
An initial condition $x(0)$. You can select this with the slider.

You will also need

A stepsize $h$. So long as this is “small enough,” the specifics don’t really matter.

How Euler works The first row of the table shows the situation at $t=0$. At that time, the value of $x$, that is $x(t=0)$ is the initial condition that you set with the slider.

In the following, whenever we write $x(t)$ we mean $x$ at the time in the last row of the table.

Knowing the value of $x(t)$ the instantaneous value of $\partial_t x$ can be found by plugging $x(t)$ into the dynamics() function.
Now that we know $\partial_t x$, we know how fast $x$ is changing. Multiply this rate of change by $h$ to get the total change of $x$ for the next step.
Add a new row to the table at $t+h$ with the $x$-value from the previous row added to the total change of $x$ from that previous row. Loop back to (a) each time the “step” button is pressed.

Select $\partial_t x = -0.5 x$ as the differential equation to solve. Press “step” several times. After each step, try to understand from the table and graphs why the new row added to the table is what it is.

Part A For $\partial_t x = -0.5 x$, which of these best describes the shape of the solution? (You will get a better picture if you set x(0) to, say, 8.)

linear decay to zero
linear growth from zero
exponential decay to zero
exponential growth from zero
exponential decay to $x = 5$
exponential growth from $x = 5$

Part B For the differential equation $\partial_t x = +0.5 x$, which of these best describes the shape of the solution? (You will get a better picture if you set x(0) to, say, 1.)

linear decay to zero
linear growth from zero
exponential decay to zero
exponential growth from zero
exponential decay to $x = 5$
exponential growth from $x = 5$

Part C For the differential equation $\partial_t x = -0.4\,(x - 5)$, which of these best describes the shape of the solution when the initial condition is $x=1$?

linear decay to zero
linear growth from zero
exponential decay to zero
exponential growth from zero
exponential decay to $x = 5$
exponential growth from $x = 5$

Part D For the differential equation $\partial_t x = -0.4\,(x - 5)$, which of these best describes the shape of the solution when the initial condition is $x=9$?

linear decay to zero
linear growth from zero
exponential decay to zero
exponential growth from zero
exponential decay to $x = 5$
exponential growth from $x = 5$

Part E For the differential equation $\partial_t x = 2\,x\,(1-x/8)$, which of these best describes the shape of the solution when the initial condition is $x=1$?

linear decay to $x=8$
exponential decay to $x=8$
exponential growth from zero followed by exponential decay to $x=8$
exponential decay to zero followed by exponential growth to $x=8$

Exercise 40.13

Figure 40.6 shows the difference between the symbolic solution to $\partial_t x = x (1-x)$ and the Euler solution. The figure shows that the Euler solution (with $dt=0.1$) has an approximation error that is small. The worst case is when $t \approx 3$ at which point it is less than 1.5 parts per million. Another good way to quantify this is the decimal place at which the calculated solutions differ.

Euler <- integrateODE(dx ~ x*(1-x), x=0.01, 
                      bounds(t=0:20), dt=0.1)
## Solution containing functions x(t).
Symb_soln = makeFun(A*exp(t)/(1 + A*exp(t)) ~ t, A=1/99)
Symb_soln(3)
## [1] 0.1686648
Euler$x(3)
## [1] 0.1686648

The approximation error occurs after the seventh decimal place.

A. How large can $dt$ be to keep the approximation error in the fourth decimal place of the answer.

B. How large can $dt$ be to keep the approximation error in the second decimal place of the answer.