Chap 40 Exercises

\[ \newcommand{\dnorm}{\text{dnorm}} \newcommand{\pnorm}{\text{pnorm}} \newcommand{\recip}{\text{recip}} \]

Exercise 1 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.

  1. \(\partial_t u = f(t)\)
  2. \(\partial_x y = g(x)\)
  3. \(\partial_t z = h(z)\)
  4. \(\partial_x y = f_2(y)\)
  5. \(\partial_t t = 1\)
  6. \(\partial_t u = g_2(u) + g_3(u)\)

Exercise 2 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)\).

  1. Write down the framework for a system of differential equations with state variables \(u, v, w\).

  2. Write down the framework for a system of differential equations with state \(z, u\).

  3. Write down the framework for a system of differential equations with state variables \(v, x, z\).

Exercise 3 Figure fig-fox-rabbit-ts 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 bump function, the rabbit graph show a little uptick near \(t=5\).

Using the R/mosaic commands given in the text to make Figure fig-fox-rabbit-ts, 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 4 We will explore some common pitfalls for using the integrateODE() and traj_plot() functions.

What is the purpose of the `integrateODE()` function?

```{mcq}
#| label: mbp204-1
#| show_hints: true
1. To translate a dynamical system into a set of numerical solutions, one for each state variable. [ correct hint: Yes  ]
2. To draw a flow field of the dynamical system. 
3. To anti-differentiate a differential equation, giving a symbolic solution. 
4. To plot the trajectory of a dynamical system. 
```
What is the proper format to use in specifying the dynamical system to `integrateODE()`?

```{mcq}
#| label: mbp204-2
#| show_hints: true
1. A set of tilde expressions, each one an individual argument. [ correct hint: Excellent!  ]
2. A set of named arguments, with the formula on the right-hand side. [ hint: You mean something like `dx = x*(1-x)`? Close, but not right. Replace the `=` by `~`. A **tilde expression** is the appropriate form. ]
3. A set of equations separated by semi-colons. [ hint: In R, arguments are always separated by **commas**. ]
4. A set of tilde expressions separated by semi-colons. [ hint: In R, arguments are always separated by **commas**. ]
```
What should the left-hand side of the tilde expression look like for a state variable `z`?

```{mcq}
#| label: mbp204-3
#| show_hints: true
1. `dz` [ correct hint: Nice.  ]
2. `z` 
3. `zdot` 
4. `dt_z` [ hint: This makes sense, and mimics the mathematical notation $\partial_t z$, but the authors of the software decided to use a shorter form. ]
```
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()`.)

```{mcq}
#| label: mbp204-4
#| show_hints: true
1. `integrateODE(dx ~ y, dy ~ -x, ...)` [ correct hint: Correct  ]
2. `integrateODE(dx = y, dy = -x, ...)` [ hint: That's  using named arguments, not the tilde expression form. ]
3. `integrateODE(dx ~ y, dy = -x, ...)` [ hint: You've got a tilde expression (correctly) for the first argument, but the second is a named argument rather than a tilde expression. It is easy to make such a mistake. ]
```
Aside from the tilde-expressions encoding the dynamical system, what other information is **absolutely essential** to specify when constructing a numerical solution?

```{mcq}
#| label: mbp204-5
#| show_hints: true
1. The initial condition for each and every state variable. [ correct hint: right-o There needs to be a starting point for the trajectory. That choice is up to you.+ ]
2. The trajectory. [ hint: This is what `integrateODE()` will calculate so long as you use it properly. ]
3. The time series. [ hint: This is what `integrateODE()` will calculate so long as you use it properly. ]
4. The initial condition for at least one of the state variables. [ hint: You need to give an initial condition for each and every state variable. ]
```
How should you specify how long you want `integrateODE()` to carry forward the solution in time, say for 10 time units?

```{mcq}
#| label: mbp204-6
#| show_hints: true
1. `bounds(t=0:10)` [ correct hint: Right  ]
2. `duration=10` 
3. `t=10` 
4. `for=10` 
```

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

  1. What are the state variables and what are the parameters?

  2. What is the dimension of the state space?

  3. 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).

  1. Based on your results in (3), what are the values of \(x\), \(y\), and \(z\) at time \(t=50\).

  2. 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\).

  3. 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.

  4. 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.

  5. 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?

  6. 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)
Loading required namespace: plotly

Exercise 6 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 7 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:

  1. 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.
  2. An initial condition \(x(0)\). You can select this with the slider.

You will also need

  1. 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.

  1. Knowing the value of \(x(t)\) the instantaneous value of \(\partial_t x\) can be found by plugging \(x(t)\) into the dynamics() function.
  2. 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.
  3. 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.

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.)

```{mcq}
#| label: odeE1-a
#| show_hints: true
1. linear decay to zero 
2. linear growth from zero 
3. exponential decay to zero [ correct hint: right-o  ]
4. exponential growth from zero 
5. exponential decay to $x = 5$ 
6. exponential growth from $x = 5$ 
```
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.)

```{mcq}
#| label: odeE1-b
#| show_hints: true
1. linear decay to zero 
2. linear growth from zero 
3. exponential decay to zero 
4. exponential growth from zero [ correct hint: Good job!  ]
5. exponential decay to $x = 5$ 
6. exponential growth from $x = 5$ 
```
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$? 

```{mcq}
#| label: odeE1-c
#| show_hints: true
1. linear decay to zero 
2. linear growth from zero 
3. exponential decay to zero 
4. exponential growth from zero 
5. exponential decay to $x = 5$ [ correct hint: You're right  ]
6. exponential growth from $x = 5$ 
```
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$? 

```{mcq}
#| label: odeE1-d
#| show_hints: true
1. linear decay to zero 
2. linear growth from zero 
3. exponential decay to zero 
4. exponential growth from zero 
5. exponential decay to $x = 5$ [ correct hint: Nice.  ]
6. exponential growth from $x = 5$ 
```
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$? 

```{mcq}
#| label: odeE1-e
#| show_hints: true
1. linear decay to $x=8$ 
2. exponential decay to $x=8$ 
3. exponential growth from zero followed by exponential decay to $x=8$ [ correct hint: Good job!  ]
4. exponential decay to zero followed by exponential growth to $x=8$ 
```

Exercise 8 Figure fig-diff-Euler-separation 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.1686645

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.

No answers yet collected