```
In [1]:
```# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *

This case study is based on "Height-Age Curves for Planted Stands of Douglas Fir, with Adjustments for Density", a working paper by Flewelling, Collier, Gonyea, Marshall, and Turnblom.

It provides "site index curves", which are curves that show the expected height of the tallest tree in a stand of Douglas firs as a function of age, for a stand where the trees are the same age.

Depending on the quality of the site, the trees might grow more quickly or slowing. So each curve is identified by a "site index" that indicates the quality of the site.

I'll start with some of the data from their Table 1. Here's the sequence of ages.

```
In [2]:
```years = [2, 3, 4, 5, 6, 8, 10, 15, 20, 25, 30,
35, 40, 45, 50, 55, 60, 65, 70]

```
Out[2]:
```

```
In [3]:
```site45 = TimeSeries([1.4, 1.49, 1.75, 2.18, 2.78, 4.45, 6.74,
14.86, 25.39, 35.60, 45.00, 53.65, 61.60,
68.92, 75.66, 81.85, 87.56, 92.8, 97.63],
index=years)

```
Out[3]:
```

Here's the series for site index 65.

```
In [4]:
```site65 = TimeSeries([1.4, 1.56, 2.01, 2.76, 3.79, 6.64, 10.44,
23.26, 37.65, 51.66, 65.00, 77.50, 89.07,
99.66, 109.28, 117.96, 125.74, 132.68, 138.84],
index=years)

```
Out[4]:
```

And for site index 85.

```
In [5]:
```site85 = Series([1.4, 1.8, 2.71, 4.09, 5.92, 10.73, 16.81,
34.03, 51.26, 68.54, 85, 100.34, 114.33,
126.91, 138.06, 147.86, 156.39, 163.76, 170.10],
index=years)

```
Out[5]:
```

Here's what the curves look like:

```
In [6]:
```plot(site85, label='SI 85')
plot(site65, label='SI 65')
plot(site45, label='SI 45')
decorate(xlabel='Time (years)',
ylabel='Height (feet)')
savefig('figs/trees-fig01.pdf')

```
```

```
In [7]:
```data = site65;

As a starting place, let's assume that the ability of the tree to gain mass is limited by the area it exposes to sunlight, and that the growth rate (in mass) is proportional to that area. In that case we can write:

$ m_{n+1} = m_n + \alpha A$

where $m_n$ is the mass of the at time step $n$, $A$ is the area exposed to sunlight, and $\alpha$ is an unknown growth parameter.

To get from $m$ to $A$, I'll make the additional assumption that mass is proportional to height raised to an unknown power:

$ m = \beta h^D $

where $h$ is height, $\beta$ is an unknown constant of proportionality, and $D$ is the dimension that relates height and mass.

We'll start by assuming $D=3$, but we'll revisit that assumption.

Finally, we'll assume that area is proportional to height squared:

$ A = \gamma h^2$

I'll specify height in feet, and choose units for mass and area so that $\beta=1$ and $\gamma=1$.

Putting all that together, we can write a difference equation for height:

$ h_{n+1}^D = h_n^D + \alpha h_n^2 $

Now let's solve it. Here's a system object with the parameters and initial conditions.

```
In [8]:
```alpha = 7
dim = 3
t_0 = get_first_label(data)
t_end = get_last_label(data)
h_0 = get_first_value(data)
system = System(alpha=alpha, dim=dim,
h_0=h_0, t_0=t_0, t_end=t_end)

```
Out[8]:
```

```
In [9]:
```def update(height, t, system):
"""Update height based on geometric model.
height: current height in feet
t: what year it is
system: system object with model parameters
"""
area = height**2
mass = height**system.dim
mass += system.alpha * area
return mass**(1/system.dim)

Test the update function with the initial conditions.

```
In [10]:
```update(h_0, t_0, system)

```
Out[10]:
```

Here's our usual version of `run_simulation`

.

```
In [11]:
```def run_simulation(system, update_func):
"""Simulate the system using any update function.
system: System object
update_func: function that computes the population next year
returns: TimeSeries
"""
results = TimeSeries()
results[system.t_0] = system.h_0
for t in linrange(system.t_0, system.t_end):
results[t+1] = update_func(results[t], t, system)
return results

And here's how we run it.

```
In [12]:
```results = run_simulation(system, update)
results.tail()

```
Out[12]:
```

Plot the results:

```
In [13]:
```def plot_results(results, data):
plot(results, ':', label='model', color='gray')
plot(data, label='data')
decorate(xlabel='Time (years)',
ylabel='Height (feet)')
plot_results(results, data)

```
```

The result converges to a straight line.

I chose the value of `alpha`

to fit the data as well as I could, but it is clear that the data have curvature that's not captured by the model.

Here are the errors:

```
In [14]:
```errors = results - data
errors.dropna()

```
Out[14]:
```

And here's the mean absolute error.

```
In [15]:
```def mean_abs_error(results, data):
return np.mean(np.abs(results-data))
mean_abs_error(results, data)

```
Out[15]:
```

This model might explain why the height of a tree grows roughly linearly:

If area is proportional to $h^2$ and mass is proportional to $h^3$, and

Change in mass is proportional to area, and

Height grows linearly, then

Area grows in proportion to $h^2$, and

Mass grows in proportion to $h^3$.

If the goal is to explain (approximate) linear growth, we might stop there. But this model does not fit the data particularly well, and it implies that trees could keep growing forever.

So we might want to do better.

As a second attempt, let's suppose that we don't know $D$. In fact, we don't, because trees are not like simple solids; they are more like fractals, which have fractal dimension.

I would expect the fractal dimension of a tree to be between 2 and 3, so I'll guess 2.5.

```
In [16]:
```alpha = 7
dim = 2.8
params = alpha, dim

```
Out[16]:
```

`System`

object.

```
In [17]:
```def make_system(params, data):
"""Makes a System object.
params: sequence of alpha, dim
data: Series
returns: System object
"""
alpha, dim = params
t_0 = get_first_label(data)
t_end = get_last_label(data)
h_0 = get_first_value(data)
return System(alpha=alpha, dim=dim,
h_0=h_0, t_0=t_0, t_end=t_end)

Here's how we use it.

```
In [18]:
```system = make_system(params, data)

```
Out[18]:
```

```
In [19]:
```def run_and_plot(alpha, dim, data):
params = alpha, dim
system = make_system(params, data)
results = run_simulation(system, update)
plot(results, ':', color='gray')

```
In [20]:
```run_and_plot(0.145, 2, data)
run_and_plot(0.58, 2.4, data)
run_and_plot(2.8, 2.8, data)
run_and_plot(6.6, 3, data)
run_and_plot(15.5, 3.2, data)
run_and_plot(38, 3.4, data)
plot(data, label='data')
decorate(xlabel='Time (years)',
ylabel='Height (feet)')

```
```

To find the parameters that best fit the data, I'll use `leastsq`

.

We need an error function that takes parameters and returns errors:

```
In [21]:
```def error_func(params, data, update_func):
"""Runs the model and returns errors.
params: sequence of alpha, dim
data: Series
update_func: function object
returns: Series of errors
"""
print(params)
system = make_system(params, data)
results = run_simulation(system, update_func)
return (results - data).dropna()

Here's how we use it:

```
In [22]:
```errors = error_func(params, data, update)

```
Out[22]:
```

`error_func`

to `leastsq`

, which finds the parameters that minimize the squares of the errors.

```
In [23]:
```best_params, details = leastsq(error_func, params, data, update)
details

```
Out[23]:
```

Using the best parameters we found, we can run the model and plot the results.

```
In [24]:
```system = make_system(best_params, data)
results = run_simulation(system, update)
plot_results(results, data)

```
```

```
In [25]:
```mean_abs_error(results, data)

```
Out[25]:
```

And the estimated fractal dimension is 3.11, which doesn't seem likely.

Let's try one more thing.

Models 1 and 2 imply that trees can grow forever, but we know that's not true. As trees get taller, it gets harder for them to move water and nutrients against the force of gravity, and their growth slows.

We can model this effect by adding a term to the model similar to what we saw in the logistic model of population growth. Instead of assuming:

$ m_{n+1} = m_n + \alpha A $

Let's assume

$ m_{n+1} = m_n + \alpha A (1 - h / K) $

where $K$ is similar to the carrying capacity of the logistic model. As $h$ approaches $K$, the factor $(1 - h/K)$ goes to 0, causing growth to level off.

Here's what the implementation of this model looks like:

```
In [26]:
```alpha = 2.0
dim = 2.5
K = 150
params = [alpha, dim, K]

```
Out[26]:
```

Here's an updated version of `make_system`

```
In [27]:
```def make_system(params, data):
"""Makes a System object.
params: sequence of alpha, dim, K
data: Series
returns: System object
"""
alpha, dim, K = params
t_0 = get_first_label(data)
t_end = get_last_label(data)
h_0 = get_first_value(data)
return System(alpha=alpha, dim=dim, K=K,
h_0=h_0, t_0=t_0, t_end=t_end)

Here's the new `System`

object.

```
In [28]:
```system = make_system(params, data)

```
Out[28]:
```

And here's the new update function.

```
In [29]:
```def update3(height, t, system):
"""Update height based on geometric model with growth limiting term.
height: current height in feet
t: what year it is
system: system object with model parameters
"""
area = height**2
mass = height**system.dim
mass += system.alpha * area * (1 - height/system.K)
return mass**(1/system.dim)

As always, we'll test the update function with the initial conditions.

```
In [30]:
```update3(h_0, t_0, system)

```
Out[30]:
```

Now we can test the error function with the new update function.

```
In [31]:
```error_func(params, data, update3)

```
Out[31]:
```

And search for the best parameters.

```
In [32]:
```best_params, details = leastsq(error_func, params, data, update3)
details

```
Out[32]:
```

With these parameters, we can fit the data much better.

```
In [33]:
```system = make_system(best_params, data)
results = run_simulation(system, update3)
plot_results(results, data)

```
```

And the mean absolute error is substantually smaller.

```
In [34]:
```mean_abs_error(results, data)

```
Out[34]:
```

The estimated fractal dimension is about 2.6, which is plausible for a tree.

Basically, it suggests that if you double the height of the tree, the mass grows by a factor of $2^{2.6}$

```
In [35]:
```2**2.6

```
Out[35]:
```

In other words, the mass of the tree scales faster than area, but not as fast as it would for a solid 3-D object.

What is this model good for?

1) It offers a possible explanation for the shape of tree growth curves.

2) It provides a way to estimate the fractal dimension of a tree based on a growth curve (probably with different values for different species).

3) It might provide a way to predict future growth of a tree, based on measurements of past growth. As with the logistic population model, this would probably only work if we have observed the part of the curve where the growth rate starts to decline.

With some help from my colleague, John Geddes, we can do some analysis.

Starting with the difference equation in terms of mass:

$m_{n+1} = m_n + \alpha A (1 - h / K) $

We can write the corresponding differential equation:

(1) $ \frac{dm}{dt} = \alpha A (1 - h / K) $

With

(2) $A = h^2$

and

(3) $m = h^D$

Taking the derivative of the last equation yields

(4) $\frac{dm}{dt} = D h^{D-1} \frac{dh}{dt}$

Combining (1), (2), and (4), we can write a differential equation for $h$:

(5) $\frac{dh}{dt} = \frac{\alpha}{D} h^{3-D} (1 - h/K)$

Now let's consider two cases:

- With infinite $K$, the factor $(1 - h/K)$ approaches 1, so we have Model 2.
- With finite $K$, we have Model 3.

Within Model 2, we'll consider two special cases, with $D=2$ and $D=3$.

With $D=2$, we have

$\frac{dh}{dt} = \frac{\alpha}{2} h$

which yields exponential growth with parameter $\alpha/2$.

With $D=3$, we have Model 1, with this equation:

$\frac{dh}{dt} = \frac{\alpha}{3}$

which yields linear growth with parameter $\alpha/3$.

This result explains why Model 1 is linear.

Within Model 3, we'll consider two special cases, with $D=2$ and $D=3$.

With $D=2$, we have

$\frac{dh}{dt} = \frac{\alpha}{2} h (1 - h/K)$

which yields logisitic growth with parameters $r = \alpha/2$ and $K$.

With $D=3$, we have

$\frac{dh}{dt} = \frac{\alpha}{3} (1 - h/K)$

which yields a first order step response; that is, it converges to $K$ like a negative exponential:

$ h(t) = c \exp(-\frac{\alpha}{3K} t) + K $

where $c$ is a constant that depends on the initial conditions.

```
In [36]:
```alpha = 10
D = 3
K = 200
params = alpha, D, K
system = make_system(params, data)
results = run_simulation(system, update3);

```
In [37]:
```t = results.index
a = alpha/3
h = (-220) * exp(-a * t / K) + K
plot(t, h)
plot(results)
decorate(xlabel='Time (years)',
ylabel='Height (feet)')

```
```

Additional resources:

Garcia, A stochastic differential equation model for the height growth of forest stands

```
In [ ]:
```