```
```{raw} latex
%%start demo%%
```
We encode the predator–prey equations via a function.
```{code-cell}
function predprey(u,p,t)
α,β = p # rename parameters for convenience
y,z = u # rename solution components
s = (y*z) / (1+β*y) # appears in both equations
return [ y*(1-α*y) - s, -z + s ]
end;
```
As before, the ODE function must accept three inputs, `u`, `p`, and `t`, even though in this case there is no explicit dependence on `t`. The second input is used to pass parameters that don't change throughout a single instance of the problem.
To specify the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable.
```{code-cell}
u₀ = [1,0.01]
tspan = (0.,60.)
α,β = 0.1,0.25
ivp = ODEProblem(predprey,u₀,tspan,[α,β])
```
You can use any `DifferentialEquations` solver on the IVP system.
```{code-cell}
sol = solve(ivp,Tsit5());
plot(sol,label=["prey" "predator"],title="Predator-prey solution")
```
We can find the discrete values used to compute the interpolated solution. The `sol.u` value is a vector of vectors.
```{code-cell}
t,u = sol.t,sol.u # extract times and solution values
@show size(u);
@show t[20];
@show u[20];
```
We can also use {numref}`Function {number} ` to find the solution.
```{code-cell}
t,u = FNC.euler(ivp,1200);
```
The solution `u` is a vector of [prey,predator] 2-vectors for each of the discrete times in `t`. Manipulating the vector-of-vectors output can be a little tricky. Here, we convert it to an $n\times 2$ matrix. Each column is one component, while each row is a single value of $t$.
```{code-cell}
u = [ u[j] for u in u, j in 1:2 ]
plot!(t[1:3:end],u[1:3:end,:],l=(1,:black),m=2,
label=["Euler prey" "Euler predator"])
```
Notice above that the accuracy of the Euler solution deteriorates rapidly.
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
When there are just two components, it's common to plot the solution in the _phase plane_, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
You can use `vars` in the plot of a solution produced by `solve` to specify the components of the solution that appear on each axis.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
plot(sol,vars=(1,2),title="Predator-prey in the phase plane",
xlabel=L"y",ylabel=L"z")
```
From this plot we can deduce that the solution approaches a periodic one, which in the phase plane is represented by a closed loop.
```{raw} html

```
```{raw} latex
%%end demo%%
```
In the rest of this chapter we present methods as though they are for scalar equations, but their application to systems is taken for granted. The generalization of error analysis can be more complicated, but our statements about order of accuracy and other properties are true for systems as well as scalars. The codes are all written to accept systems.
## Transformation of high-order systems
Fortunately, the ability to solve first-order ODE systems implies the ability to solve systems of higher differential order, too. The reason is that there is a systematic way to turn a higher-order problem into a first-order one of higher dimension.
(example-nlosc3)=
````{proof:example}
Consider the nonlinear initial-value problem
```{math}
y''+(1+y')^3 y = 0, \qquad y(0)= y_0, \quad y'(0) = 0.
```
In order to write this problem as a first-order system we define two scalar unknown functions, $u_1 = y$ and $u_2 = y'$. With these definitions, we have the two differential equations
```{math}
\begin{split}
u_1' &= u_2, \\
u_2' &= -(1+u_2)^3 u_1,
\end{split}
```
which is a first-order system in two dimensions. The initial
condition of the system is
```{math}
u_1(0) = y_0, \quad u_2(0) = 0.
```
````
(example-systems-coupledpendula)=
````{proof:example}
Two identical pendulums suspended from the same rod and swinging in parallel planes can be modeled as the second-order system
```{math}
\begin{split}
\theta_1''(t) +\gamma \theta_1' + \frac{g}{L} \sin \theta_1 +
k(\theta_1-\theta_2) &= 0,\\
\theta_2''(t) +\gamma \theta_2' + \frac{g}{L} \sin \theta_2 +
k(\theta_2-\theta_1) &= 0,
\end{split}
```
where $\theta_1$ and $\theta_2$ are angles made by the two pendulums, $L$ is the length of each pendulum, $\gamma$ is a frictional parameter, and $k$ is a parameter describing a torque produced by the rod when it is twisted. We can convert this problem into a first-order system using the substitutions
```{math}
u_1 = \theta_1, \quad u_2 = \theta_2, \quad u_3 = \theta_1', \quad
u_4 = \theta_2'.
```
With these definitions the system becomes
```{math}
\begin{split}
u_1' &= u_3, \\
u_2' &= u_4, \\
u_3' &= -\gamma u_3 - \frac{g}{L}\sin u_1 + k(u_2-u_1), \\
u_4' &= -\gamma u_4 - \frac{g}{L}\sin u_2 + k(u_1-u_2),
\end{split}
```
which is a first-order system in four dimensions. To complete the description of the problem, you need to specify values for $\theta_1(0)$, $\theta_1'(0)$, $\theta_2(0)$, and $\theta_2'(0)$.
````
The trick illustrated in the preceding examples is always available. Suppose $y$ is a scalar dependent variable in the system. You should introduce a component of $\mathbf{u}$ for $y$, $y'$, etc., up to but not including the highest derivative appearing anywhere for $y$. This is done for each scalar variable in the original system. There should be one component of $\mathbf{u}$ for each scalar initial condition given. Many equations for the first-order system then come from the trivial relationships among all the lower derivatives. The remaining equations for the system come from the original, high-order equations. In the end, there must be as many scalar component equations as unknown first-order variables.
(demo-systems-coupledpendula)=
```{proof:demo}
```
```{raw} html
```
```{raw} latex
%%start demo%%
```
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
Let's implement the coupled pendulums from {numref}`Example {number} `. The pendulums will be pulled in opposite directions and then released together from rest.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
The `similar` function creates an array of the same size and type as a given value, without initializing the contents.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
function couple(u,p,t)
γ,L,k = p
g = 9.8
udot = similar(u)
udot[1:2] .= u[3:4]
udot[3] = - γ*u[3] - (g/L)*sin(u[1]) + k*(u[2]-u[1])
udot[4] = - γ*u[4] - (g/L)*sin(u[2]) + k*(u[1]-u[2])
return udot
end
u₀ = [1.25,-0.5,0,0]
tspan = (0.,50.);
```
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
First we check the behavior of the system when the pendulums are uncoupled, i.e., when $k=0$.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
Here `vars` is used to plot two components as functions of time.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
γ,L,k = 0,0.5,0
ivp = ODEProblem(couple,u₀,tspan,[γ,L,k])
sol = solve(ivp,Tsit5())
plot(sol,vars=[1,2],label=[L"\theta_1" L"\theta_2"],
xlims=[20,50],title="Uncoupled pendulums")
```
You can see that the pendulums swing independently. Because the model is nonlinear and the initial angles are not small, they have slightly different periods of oscillation, and they go in and out of phase.
With coupling activated, a different behavior is seen.
```{code-cell}
k = 1
ivp = ODEProblem(couple,u₀,tspan,[γ,L,k])
sol = solve(ivp,Tsit5())
plot(sol,vars=[1,2],label=[L"\theta_1" L"\theta_2"],
xlims=[20,50],title="Coupled pendulums")
```
The coupling makes the pendulums swap energy back and forth.
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Exercises
1. ✍ Rewrite the given higher order problems as first-order systems.
**(a)** $y'''-3y''+3 y' -y = t, \: y(0) = 1, \: y'(0) = 2, \: y''(0) = 3$
**(b)** $y'' + 4 (x^2-1)y' + y = 0, \: y(0) = 2, \: y'(0) = -1$
**(c)** For a given constant $a$,
```{math}
\begin{split}
x'' + \frac{a x}{(x^2+y^2)^{3/2}} &= 0,\\
y'' + \frac{a y}{(x^2+y^2)^{3/2}} &= 0,
\end{split}
```
with initial values $x(0) = 1$, $x'(0)=y(0) = 0$, $y'(0)=3$
**(d)** $y^{(4)} -y = e^{-t}, \: y(0) = 0, \: y'(0) = 0, \: y''(0) = 1,\: y'''(0) = 0$
**(e)** $y'''-y''+y'-y = t, \: y(0) = 1, \: y'(0) = 2, \: y''(0) = 3$
2. ✍ Write the given IVP as a system. Then do two steps of Euler's method by hand (perhaps with a calculator) with the indicated step size $h$. Using the given exact solution, compute the error after the second step.
**(a)** $y''+ 4y = 4t, \: y(0) = 1,\: y'(0) = 1; \: \hat{y}(t) = t+\cos (2t),\: h=0.1$
**(b)** $y''- 4y = 4t, \: y(0) = 2,\: y'(0) = -1; \: \hat{y}(t) = e^{2t} + e^{-2t}-t,\: h=0.1$
**(c)** $2 x^2 y'' +3xy' - y = 0, \: y(2) = 1, \: y'(2) = -1/2, \: \hat{y}(x) = 2/x, h = 1/8$
**(d)** $2 x^2 y'' +3xy' - y = 0,\: y(1) = 4, \: y'(1) = -1, \: \hat{y}(x) = 2(x^{1/2} + x^{-1}), h=1/4$
3. ⌨ Solve the following IVPs using {numref}`Function {number}