Homework 6
1. Given n scalar measurements {x1 , x2 , . . . , xN }, you can estimate the mean using the formula 1 N
x ̄=N xi i=1
and the standard deviation using the formula:
1 N
( x i − x ̄ ) 2
One of the more puzzling aspects of statistics, at least to me, is that you normalize by N − 1 rather than by N when estimating the standard deviation. The justification for this is actually quite complicated and often a simpler explanation is given regarding the number of degrees of freedom (e.g. you lose one when you estimate the mean).
Using Python, show that the formula for s gives a better estimate of the standard deviation
s = N − 1
i=1
than the formula:
1 N
s 1 = N
( x i − x ̄ ) 2
i=1
To do this, estimate the standard deviation for samples of size N = 2…40. Assume that you are sampling from a normal distribution with a mean of 2 and standard deviation of 5. To get a good estimate for the standard deviation, repeat this computer experiment 50,000 times and then plot the average for the estimate of the standard deviation using the different formulas.
In addition, show that following formula does an even better job of estimating the standard deviation, particularly when N is small:
1 N
s 2 = N − 1 . 5
As background, we use the notation x ∼ N (0, 1) to say that x is normally distributed with a mean of 0 and a standard deviation of 1. If y = ax+b and x ∼ N(0,1), then y ∼ N(a,b). In other words, you can generate a normally distributed random variable with mean a and standard deviation b in Python with the command: b*np.random.randn()+a.
1
( x i − x ̄ ) 2
i=1
Your results should look like:
Homework 6
5.5 5.0 4.5 4.0 3.5 3.0
N-1
N
N-1.5
Exact
5 10 15 20 25 30 35 40 N
2. The following differential equations, known as the Lokta-Volterra equations, are used to model predator-prey dynamics:
dx = αx − βxy, dt
dy =δxy−γy, dt
where x denotes the number of prey and y the number of predators (these values are dimen- sionless). The variables α, β, δ, and γ are parameters. Use Python to solve the Lokta-Volterra equations for the parameter values α = β = δ = γ = 1. Use the following initial conditions: x(0) = 1 and y(0) = 0.1. Your results should look like
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
0 5 10 15 20 25 time
Prey
Predator
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
0.0 0.5
1.0 1.5 2.0 prey
2.5 3.0 3.5
The second plot is what is known as a phase portrait, where y(t) is plotted against x(t). 3. The following differential equation is used to model the oscillatory motion of a pendulum:
d2θ + g sin(θ) = 0, dt2 l
2
predator
STD
Homework 6
where g is the gravitational constrant (9.8 m2/s), l is the length of the pendulum in me- ters, and θ the angular displacement. See https://en.wikipedia.org/wiki/Pendulum_ (mathematics) for further information.
You task is to determine how the period of these oscillations varies as the length of the pendulum. The initial conditions should be: θ(0) = 90 degrees and θ ̇(0) = 0 degrees/s. The results should look like:
20 15 10
5
0 20 40 60 80 100 length (m)
4. Consider a batch reactor where the reversible reaction occurs:
r1
A −→ B, r1 = 25200 exp(−41868/RT )[A] ([=] M/hr)
r2 7
B −→ A, r2 = 1.8 × 10 exp(−83737/RT )[B] ([=] M/hr)
where R = 8.314. Determine the temperature where the concentration of B is greatest after 24 hrs? Plot the concentration profile of B at this optimal temperature. Assume the initial concentration is: [A](0) = 1 and [B](0) = 0. The solution is plotted below. The max temperature is approximately 450◦K (the answer will depend on how fine of grid you use). The right figure shows the concentration profiles at the optimal temperature.
1.0 0.8 0.6 0.4 0.2 0.0
0 200
400 600 Temperature
800 1000
1.0 0.8 0.6 0.4 0.2 0.0
0 5 10 15 20 25 time
A
B
3
B
concentration
period (s)
5. Solve the Lorenz equations
Homework 6
dx = σ(y − x) dt
dy =rx−y−xz dt
dz =xy−bz dt
for the parameter values: σ = 10, b = 8/3, and r = 28. The initial conditions should be: x(0) = 1, y(0) = 0, z(0) = 0. Solve the equations over the time interval: [0,100]. These equations provide an example of chaotic dynamics. Plot the solution for x(t) versus y(t). The plot should look like:
4