Homework # 3
ME526/NSE526
Due: November 5
1. You are encouraged to work in a group of up to 2 students and submit one solution per
group.
2. Your solution must be clearly legible. Illegible work may not be graded and returned without
any points. Although not necessary, you may type your work.
3. All problems must be solved. However, all problems may not be graded. A random sample
of problems will be selected for grading.
4. If you are required to write a computer program, attach your code with several comment
statements on the code wherever possible.
1
1. The following scheme is proposed to solve y′ = f(y):
y∗ = yn + γ1hf(yn) (1)
y∗∗ = y∗ + γ2hf(y∗) + ω2hf(yn) (2)
yn+1 = y
∗∗ + γ3hf(y∗∗) + ω3hf(y∗) (3)
where γ1 = 8/15, γ2 = 5/12, γ3 = 3/4, ω2 = −17/60, and ω3 = −5/12 with h being the time
step.
(a) What is the order of accuracy of this method? Prove the order by applying the method
to y′ = λy and obtaining a finite difference approximation that can be written in the
yn+1 = σyn.
(b) Draw a stability diagram in the (hλR, hλI) plane for this method applied to the model
problem y′ = λy. What is the maximum step size for λ pure imaginary and λ negative
real?
2. Consider motion of a simple pendulum (see Fig 1) consisting of a mass m attached to a
string of length ` and immersed in a viscous fluid adding some damping resistance. The
governing equation is given by
θ′′ + cθ′ +
g
`
sin(θ) = 0, (4)
where g is the gravitational acceleration, and θ the angle made by the string with the vertical
axis, c is the damping coefficient. For small angles θ, the linearized equation of motion is:
θ′′ + cθ′ +
g
`
θ = 0. (5)
Let g = 9.81 m/sec2, and ` = 0.6 m, c = 4 sec−1. Assume that the pendulum starts from
rest (i.e. θ′(t = 0) = 0) with θ(t = 0) = 100.
Figure 1: Schematic of a simple pendulum of mass m and length ` immersed in a viscous fluid.
(a) Evaluate the exact solution for the linearized system (equation 5).
2
(b) Consider the following numerical schemes to advance the solution over the time interval
0 ≤ t ≤ 6. For each scheme clearly indicate the system of discretized equations (to
advance from t = tn to t = tn+1) along with the proper initial conditions to be used.
i. The Forward Euler,
ii. The Second Order Runge-Kutta,
iii. The Fourth Order Runge-Kutta.
(c) Perform linear stability analysis and discuss quantitatively and in detail the stability
limits for each scheme.
(d) Evaluate the solution to the linearized problem 1, θ(t), using all the schemes for 0 ≤
t ≤ 6. Use ∆t = 0.15, 0.5, 1 s. For each case, and on separate figures, plot the solution
over the time period together with the exact solution to the linearized system2 Use
subplot to plot all figures (that is total of 9 figures, 3 for each scheme and 3
different time-steps) on a single page. Clearly label all plots.
(e) Consider the non-linear equation 4. For an undamped system (i.e. c = 0), solve the
differential equations with θ(t = 0) = 600 with RK4. What steps have you taken to be
certain of the accuracy of your results? Why should your results be believable? Plot
the solution for 0 ≤ t ≤ 6. Compare your solution with exact solution of the linearized
undamped governing equation (i.e. solution to equation 5 with c = 0). Comment on
what is the effect of the non-linearity on the solution. For the numerical approach, find
out any limit on the time-step to obtain stable solution for the non-linear problem using
RK4 (This should be done by increasing the time-steps in small increments and finding
a time-step beyond which unphysical solution is obtained. List down the time-step at
which this happens). Compare this time-step obtained from numerical experiments to
the time-step restriction provided by the linear-stability theory for the undamped
problem. How do they compare? Comment on your findings.
3. Consider chemical reactions during food digestion in our bodies. An enzyme (E) combines
with a substrate (S) to form a complex (ES). The complex can dissociate back to E and S
or it can proceed to form product P. The time history of the reactions
E + S
k1−⇀↽−
k2
ES
k3→ E + P, (6)
is gorverned by the following rate equations The above reaction are governed by the following
1You will have to write a computer program to do this with separate routines for each scheme.
2In your plot θ should be in degrees. Make sure you use separate line-styles or symbols to distinguish the plots,
also clearly label x and y axes, and use the same plotting format for all figures (e.g. the axis range etc.). If your
answer blows up, plot the answer up to a certain range of θ axis to allow comparison with other schemes.
3
rate equations:
dCS
dt
= −k1CSCE + k2CES (7)
dCE
dt
= −k1CSCE + (k2 + k3)CES (8)
dCES
dt
= k1CSCE − (k2 + k3)CES (9)
dCP
dt
= k3CES (10)
where Cis are concentrations and k1, k2 and k3 are reaction rate constants given as
k1 = 2× 103; k2 = 1× 10−3; k3 = 10.0. (11)
Initially, CS(0) ∼ 1.0, CE(0) = 5×10−5, CES(0) = 0, and CP = 0. Using C = (C1, C2, C3, C4) =
(CS, CE, CES, CP), Solve the above system numerically to steady state (let t = 2500s repre-
sent steady state) using the following schemes.
(a) RK4. Do you think ∆t = 0.0005 will give you a stable solution? Justify your answer
based on Stability Analysis at t = 0. Use this time-step to obtain the steady state
solution. Note the computing time required to obtain steady state. Compare steady
state solution to the exact. Plot all concentrations on the same plot (use axes ranges
from 10−10 to 1.5× 100 for concentrations and 10−5 to 3000 for time. Provide solutions
in a log-log plot as well as semi-log plot (where the concentrations are on a log scale
and time in a linear scale).
(b) Setup the problem with a linearized trapezoidal scheme. To do this, first write down
the finite-difference equations for the standard trapezoidal scheme. Since this is an
implicit method, you will have non-linear terms for the unknown concentrations at
time level tn+1. Apply the linearization technique to the unknown non-linear terms in
your finite-difference approximation (fda) and re-write the fda indicating the time level
at which each term is evaluated. Your final answer should be in the form: MY = NX+b
where M, N are matrices, Y , X, and b are some vectors. Clearly indicate the elements
of each matrix and vector and also the time level at which each element is evaluated.
Would a time step of 0.001 provide an accurate solution at steady state? Solve and
check. What advantages would such a linearized central differencing scheme would have
over say the fourth-order RK method?
(c) Compare the computing time required for the three approaches above.
4
