辅导案例-ME526/NSE526

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. 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θ′ +

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θ′ +

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

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

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−⇀↽−

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.

rate equations:

dCS

= −k1CSCE + k2CES (7)

dCE

= −k1CSCE + (k2 + k3)CES (8)

dCES

= k1CSCE − (k2 + k3)CES (9)

dCP

= 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?