Exercises

27.9 EXERCISES
27.9 Exercises
27.1
27.2
27.3
Use an iteration procedure to find the root of the equation 40x = exp x to four
significant figures.
Using the Newton–Raphson procedure find, correct to three decimal places, the
root nearest to 7 of the equation 4x3 + 2x2 − 200x − 50 = 0.
Show the following results about rearrangement schemes for polynomial equations.
(a) That if a polynomial equation g(x) ≡ xm − f(x) = 0, where f(x) is a
polynomial of degree less than m and for which f(0) = 0, is solved using
a rearrangement iteration scheme xn+1 = [ f(xn )]1/m , then, in general, the
scheme will have only first-order convergence.
(b) By considering the cubic equation
x3 − ax2 + 2abx − (b3 + ab2 ) = 0
for arbitrary non-zero values of a and b, demonstrate that, in special cases, the
same rearrangement scheme can give second- (or higher-) order convergence.
27.4
27.5
The square root of a number N is to be determined by means of the iteration
scheme
xn+1 = xn 1 − N − x2n f(N) .
Determine how
√ the process has second-order convergence.
√to choose f(N) so that
Given that 7 ≈ 2.65, calculate 7 as accurately as a single application of the
formula will allow.
Solve the following set of simultaneous equations using Gaussian elimination
(including interchange where it is formally desirable):
x1 + 3x2 + 4x3 + 2x4
2x1 + 10x2 − 5x3 + x4
4x2 + 3x3 + 3x4
−3x1 + 6x2 + 12x3 − 4x4
27.6
The following table of values of a polynomial p(x) of low degree contains an
error. Identify and correct the erroneous value and extend the table up to x = 1.2.
x
0.0
0.1
0.2
0.3
0.4
27.7
= 0,
= 6,
= 20,
= 16.
p(x)
0.000
0.011
0.040
0.081
0.128
x
0.5
0.6
0.7
0.8
0.9
p(x)
0.165
0.216
0.245
0.256
0.243
Simultaneous linear equations that result in tridiagonal matrices can sometimes
be treated as three-term recurrence relations, and their solution may be found
in a similar manner to that described in chapter 15. Consider the tridiagonal
simultaneous equations
xi−1 + 4xi + xi+1 = 3(δi+1,0 − δi−1,0 ),
i = 0, ±1, ±2, . . . .
Prove that, for i > 0, the equations have a general solution of the form xi =
αpi + βq i , where p and q are the roots of a certain quadratic equation. Show that
a similar result holds for i < 0. In each case express x0 in terms of the arbitrary
constants α, β, . . . .
Now impose the condition that xi is bounded as i → ±∞ and obtain a unique
solution.
1033
NUMERICAL METHODS
27.8
A possible rule for obtaining an approximation to an integral is the mid-point
rule, given by
x0 +∆x
f(x) dx = ∆x f(x0 + 12 ∆x) + O(∆x3 ).
x0
27.9
27.10
Writing h for ∆x, and evaluating all derivates at the mid-point of the interval
(x, x + ∆x), use a Taylor series expansion to find, up to O(h5 ), the coefficients of
the higher-order errors in both the trapezium and mid-point rules. Hence find a
linear combination of these two rules that gives O(h5 ) accuracy for each step ∆x.
Although it can easily be shown, by direct calculation, that
∞
1
e−x cos(kx) dx =
,
1 + k2
0
the form of the integrand is appropriate for Gauss–Laguerre numerical integration. Using a 5-point formula, investigate the range of values of k for which the
formula gives accurate results. At about what value of k do the results become
inaccurate at the 1% level?
Using the points and weights given in table 27.9, answer the following questions.
(a) A table of unnormalised Hermite polynomials Hn (x) has been spattered with
ink blots and gives H5 (x) as 32x5 −?x3 + 120x and H4 (x) as ?x4 −?x2 + 12,
where the coefficients marked ? cannot be read. What should they read?
(b) What is the value of the integral
∞
2
e−2x
I=
dx,
2
−∞ 4x + 3x + 1
as given by a 7-point integration routine?
27.11
Consider the integrals Ip defined by
Ip =
1
x2p
√
dx.
1 − x2
−1
(a) By setting x = sin θ and using the results given in exercise 2.42, show that Ip
has the value
2p − 1 2p − 3
1 π
Ip = 2
···
.
2p 2p − 2
2 2
(b) Evaluate Ip for p = 1, 2, . . . , 6 using 5- and 6-point Gauss–Chebyshev integration (conveniently run on a spreadsheet such as Excel) and compare the
results with those in (a). In particular, show that, as expected, the 5-point
scheme first fails to be accurate when the order of the polynomial numerator
(2p) exceeds (2×5)−1 = 9. Likewise, verify that the 6-point scheme evaluates
I5 accurately but is in error for I6 .
27.12
In normal use only a single application of n-point Gaussian quadrature is made,
using a value of n that is estimated from experience to be ‘safe’. However, it is
instructive to examine what happens when n is changed in a controlled way.
(a) Evaluate the integral
5
In =
√
7x − x2 − 10 dx
2
using n-point Gauss–Legendre formulae for n = 2, 3, . . . , 6. Estimate (to 4
s.f.) the value I∞ you would obtain for very large n and compare it with the
result I obtained by exact integration. Explain why the variation of In with
n is monotonically decreasing.
1034
27.9 EXERCISES
(b) Try to repeat the processes described in (a) for the integrals
5
1
√
Jn =
dx.
7x − x2 − 10
2
Why is it very difficult to estimate J∞ ?
27.13
Given a random number η uniformly distributed on (0, 1), determine the function
ξ = ξ(η) that would generate a random number ξ distributed as
(a) 2ξ√ on 0 ≤ ξ < 1,
(b) 32 ξ on 0 ≤ ξ < 1,
π
πξ
(c)
cos
on − a ≤ ξ < a,
4a
2a
(d) 12 exp(− | ξ |) on − ∞ < ξ < ∞.
27.14
27.15
27.16
27.17
A, B and C are three circles of unit radius with centres in the xy-plane at
(1, 2), (2.5, 1.5) and (2, 3), respectively. Devise a hit or miss Monte Carlo calculation
to determine the size of the area that lies outside C but inside A and B, as well
as inside the square centred on (2, 2.5), that has sides of length 2 parallel to the
coordinate axes. You should choose your sampling region so as to make the
estimation as efficient as possible. Take the random number distribution to be
uniform on (0, 1) and determine the inequalities that have to be tested using the
random numbers chosen.
Use a Taylor series to solve the equation
dy
+ xy = 0,
y(0) = 1,
dx
evaluating y(x) for x = 0.0 to 0.5 in steps of 0.1.
Consider the application of the predictor–corrector method described near the
end of subsection 27.6.3 to the equation
dy
= x + y,
y(0) = 0.
dx
Show, by comparison with a Taylor series expansion, that the expression obtained
for yi+1 in terms of xi and yi by applying the three steps indicated (without any
repeat of the last two) is correct to O(h2 ). Using steps of h = 0.1 compute the
value of y(0.3) and compare it with the value obtained by solving the equation
analytically.
A more refined form of the Adams predictor–corrector method for solving the
first-order differential equation
dy
= f(x, y)
dx
is known as the Adams–Moulton–Bashforth scheme. At any stage (say the nth)
in an Nth-order scheme, the values of x and y at the previous N solution points
are first used to predict the value of yn+1 . This approximate value of y at the
next solution point, xn+1 , denoted by ȳn+1 , is then used together with those at the
previous N − 1 solution points to make a more refined (corrected) estimation of
y(xn+1 ). The calculational procedure for a third-order scheme is summarised by
the two following two equations:
ȳn+1 = yn + h(a1 fn + a2 fn−1 + a3 fn−2 )
yn+1 = yn + h(b1 f(xn+1 , ȳn+1 ) + b2 fn + b3 fn−1 )
(predictor),
(corrector).
(a) Find Taylor series expansions for fn−1 and fn−2 in terms of the function
fn = f(xn , yn ) and its derivatives at xn .
1035
NUMERICAL METHODS
(b) Substitute them into the predictor equation and, by making that expression
for ȳn+1 coincide with the true Taylor series for yn+1 up to order h3 , establish
simultaneous equations that determine the values of a1 , a2 and a3 .
(c) Find the Taylor series for fn+1 and substitute it and that for fn−1 into the
corrector equation. Make the corrected prediction for yn+1 coincide with the
true Taylor series by choosing the weights b1 , b2 and b3 appropriately.
(d) The values of the numerical solution of the differential equation
dy
2(1 + x)y + x3/2
=
dx
2x(1 + x)
at three values of x are given in the following table:
x
y(x)
0.1
0.030 628
0.2
0.084 107
0.3
0.150 328
Use the above predictor–corrector scheme to find the value of y(0.4) and
compare your answer with the accurate value, 0.225 577.
27.18
If dy/dx = f(x, y) then show that
∂2 f
∂2 f
∂f ∂f
∂2 f
d2 f
= 2 + 2f
+ f2 2 +
+f
dx2
∂x
∂x∂y
∂y
∂x ∂y
∂f
∂y
2
.
Hence verify, by substitution and the subsequent expansion of arguments in
Taylor series of their own, that the scheme given in (27.79) coincides with the
Taylor expansion (27.68), i.e.
yi+1 = yi + hyi(1) +
27.19
h2 (2) h3 (3)
y + yi + · · · ,
2! i
3!
up to terms in h3 .
To solve the ordinary differential equation
du
= f(u, t)
dt
for f = f(t), the explicit two-step finite difference scheme
un+1 = αun + βun−1 + h(µfn + νfn−1 )
may be used. Here, in the usual notation, h is the time step, tn = nh, un = u(tn )
and fn = f(un , tn ); α, β, µ, and ν are constants.
(a) A particular scheme has α = 1, β = 0, µ = 3/2 and ν = −1/2. By considering
Taylor expansions about t = tn for both un+j and fn+j , show that this scheme
gives errors of order h3 .
(b) Find the values of α, β, µ and ν that will give the greatest accuracy.
27.20
Set up a finite difference scheme to solve the ordinary differential equation
d2 φ dφ
+
=0
dx2
dx
in the range 1 ≤ x ≤ 4, subject to the boundary conditions φ(1) = 2 and
dφ/dx = 2 at x = 4. Using N equal increments, ∆x, in x, obtain the general
difference equation and state how the boundary conditions are incorporated
into the scheme. Setting ∆x equal to the (crude) value 1, obtain the relevant
simultaneous equations and so obtain rough estimates for φ(2), φ(3) and φ(4).
Finally, solve the original equation analytically and compare your numerical
estimates with the accurate values.
x
1036
27.9 EXERCISES
27.21
Write a computer program that would solve, for a range of values of λ, the
differential equation
dy
1
,
= dx
x2 + λy 2
27.22
y(0) = 1,
using a third-order Runge–Kutta scheme. Consider the difficulties that might
arise when λ < 0.
Use the isocline approach to sketch the family of curves that satisfies the nonlinear first-order differential equation
dy
a
.
= dx
x2 + y 2
27.23
For some problems, numerical or algebraic experimentation may suggest the
form of the complete solution. Consider the problem of numerically integrating
the first-order wave equation
∂u
∂u
+A
= 0,
∂t
∂x
in which A is a positive constant. A finite difference scheme for this partial
differential equation is
u(p, n) − u(p − 1, n)
u(p, n + 1) − u(p, n)
+A
= 0,
∆t
∆x
where x = p∆x and t = n∆t, with p any integer and n a non-negative integer.
The initial values are u(0, 0) = 1 and u(p, 0) = 0 for p = 0.
(a) Carry the difference equation forward in time for two or three steps and
attempt to identify the pattern of solution. Establish the criterion for the
method to be numerically stable.
(b) Suggest a general form for u(p, n), expressing it in generator function form,
i.e. as ‘u(p, n) is the coefficient of sp in the expansion of G(n, s)’.
(c) Using your form of solution (or that given in the answers!), obtain an
explicit general expression for u(p, n) and verify it by direct substitution into
the difference equation.
(d) An analytic solution of the original PDE indicates that an initial disturbance propagates undistorted. Under what circumstances would the difference scheme reproduce that behaviour?
27.24
In exercise 27.23 the difference scheme for solving
∂u
∂u
+
= 0,
∂t
∂x
in which A has been set equal to unity, was one-sided in both space (x) and
time (t). A more accurate procedure (known as the Lax–Wendroff scheme) is
u(p + 1, n) − u(p − 1, n)
u(p, n + 1) − u(p, n)
+
∆t
2∆x
∆t u(p + 1, n) − 2u(p, n) + u(p − 1, n)
.
=
2
(∆x)2
(a) Establish the orders of accuracy of the two finite difference approximations
on the LHS of the equation.
(b) Establish the accuracy with which the expression in the brackets approximates ∂2 u/∂x2 .
(c) Show that the RHS of the equation is such as to make the whole difference
scheme accurate to second order in both space and time.
1037
NUMERICAL METHODS
27.25
Laplace’s equation,
∂2 V
∂2 V
+
= 0,
∂x2
∂y 2
is to be solved for the region and boundary conditions shown in figure 27.7.
V = 80
−∞
40
40
40
40
40
20
20
20
40
40
∞
V =0
Figure 27.7 Region, boundary values and initial guessed solution values.
27.26
Starting from the given initial guess for the potential values V , and using the
simplest possible form of relaxation, obtain a better approximation to the actual
solution. Do not aim to be more accurate than ± 0.5 units, and so terminate the
process when subsequent changes would be no greater than this.
Consider the solution, φ(x, y), of Laplace’s equation in two dimensions using a
relaxation method on a square grid with common spacing h. As in the main text,
denote φ(x0 + ih, y0 + jh) by φi,j . Further, define φm,n
i,j by
φm,n
i,j ≡
∂m+n φ
∂xm ∂y n
evaluated at (x0 + ih, y0 + jh).
(a) Show that
2,2
0,4
φ4,0
i,j + 2φi,j + φi,j = 0.
(b) Working up to terms of order h5 , find Taylor series expansions, expressed in
terms of the φm,n
i,j , for
S±,0 = φi+1,j + φi−1,j ,
S0,± = φi,j+1 + φi,j−1 .
(c) Find a corresponding expansion, to the same order of accuracy, for φi±1,j+1 +
φi±1,j−1 and hence show that
S±,± = φi+1,j+1 + φi+1,j−1 + φi−1,j+1 + φi−1,j−1
has the form
h4 4,0
0,4
(φ + 6φ2,2
i,j + φi,j ).
6 i,j
(d) Evaluate the expression 4(S±,0 + S0,± ) + S±,± and hence deduce that a possible
relaxation scheme, good to the fifth order in h, is to recalculate each φi,j as
the weighted mean of the current values of its four nearest neighbours (each
1
with weight 15 ) and its four next-nearest neighbours (each with weight 20
).
2,0
0,2
2
4φ0,0
i,j + 2h (φi,j + φi,j ) +
1038