ENGR20005
Numerical Methods in Engineering
Workshop 7
Part A: MATLAB Livescripts
7.1 The livescript ENGR20005 Workshop7p1.mlx runs through the use of MATLAB func- tions to integrate functions.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Modify the livescript to compute the integral
10 0
where f(x) is given by i. f(x)=x3 +6x+9
ii. f(x) = sin2 (x) tan (x2)
7.2 The livescript ENGR20005 Workshop7p2.mlx runs through numerical integration.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Modify the livescript to evaluate the integrals in Question 7.1.
1
f(x) dx
(7.1)
Part B: Problems
7.3 Consider the function
f(x) = 6×3 +x2 −11x+4 What is the error in evaluating the integral
1
f(x) dx
0
using both the Trapezoidal and Simpson’s rule with a single interval.
7.4 * In applied mathematics, the error function
1x
exp(−τ2)dτ
often appears in statistics as well as phenomena involving diffusion.
(7.2)
erf(x)= √π
(7.3)
(a) Write a MATLAB function (.m file) that computes Eq. (7.3), using Gaussian– Legendre quadrature with k points. Your function should take x and k as inputs, and output the value of erf(x).
To obtain the quadrature points and weights, the function gauss legendre points.m has been provided.
Compare your function with MATLAB’s inbuilt erf() function.
(b) The evolution of a flow created by a infinitely long plate accelerating instantly
from rest to a velocity of U is given by
u = U1 − erf√y (7.4)
4νt
Use your answer to (a) to plot u within the domain (y, t) ∈ [0, 2] × [0, 5], assuming
a viscosity of ν = 1 m2/s and a plate velocity of U = 1 m/s.
7.5 While Gauss–Legendre quadrature is useful most of the time, sometimes it’s necessary
to include the end nodes at x = −1, 1.
Assume that the quadrature formula is given by
1 −1
(a) Determine the quadrature weights and points for the two, three, four, and five- point quadrature formulas.
I =
f(x)dx=w0f(−1)+w1f(x1)+···+wn−1f(xn−1)+wnf(1) (7.5)
2
−x
(b) (c)
What is the highest order of polynomial that Eq. (7.5) will integrate exactly. i. Show that the Legendre differential equation may be written as
d (t2 − 1)dPk (t) = k(k + 1)Pk(t) dt dt
Pk(τ)dτ
−1 which are known as the Completed Lobatto polynomials.
iii. Determine the Complete Lobatto polynomials when k = 1, 2, 3, 4.
iv. Determine the roots of the polynomials you found in part iii.
(7.6)
ii. Integrate Eq. (7.6) once and show that
ˆ 2 dPk t
Lok+1(t)=(1−t ) dt =−k(k+1)
(d)
7.6 [Challenge.]One technique of solving linear differential equations of the form
aN f (N ) (t) + aN −1 f (N −1) (t) + · · · + a1 f (1) (t) + a0 f (t) = b(t) is to take the Laplace transform of the differential equation using
0
s-space. You’ll learn more about this in MAST20029 Engineering Mathematics.
Our goal here is to create some functions to perform Laplace transforms and their
inverses.
(a) Write a function (.m file) that applies the Trapezoidal rule to evaluate Eq. (7.8).
Use what you found in part (c) to determine higher order Gauss–Lobatto quadra- ture rules.
ˆ∞ f(s) =
f(t)e−st dt
This essentially transforms a calculus problem in t-space to an algebraic problem in
Your function should take an anonymous function and a vector of s values as an ˆ
input, and output a vector of f’s at s.
Hint: To numerically compute integrals defined over an infinite domain, we note
that for a convergent integral
(b) Use your function to evaluate the Laplace transforms of 3
∞LL f(t)dt= lim f(t)dt≈
f(t)dt
So if we choose L to be large enough such that the integrand is close to zero, the
0 L→∞0 0 approximation will be close enough for practical purposes.
(7.7)
(7.8)
i. f(t) = e−t
ii. f(t)=sinh(t)
Compare your function with MATLAB’s inbuilt laplace() function.
(c) Under the assumptions that f(t) is non-negative, f(t) = 0 for t < 0, the real part (R) of fˆ is even, and the imaginary part (I) of fˆ is odd, the inverse Laplace transform is given by
2eγt ∞ ˆ
ˆ
f, which we’ll just guess and hope the integral converges to the correct value.
Apply the Trapezoidal rule with equally spaced nodes to Eq. (7.9) and show that
∆eγtˆ 2∆eγt∞ˆ
f(t)≈ π R f(γ+iω) + π
(d) Assume that ∆ = π 2t
2t
A eA/2 ∞
f(t)≈
R fˆ +
(−1)kR fˆ
2t
R f(γ+iω) cos(ωt)dω (7.9)
f(t)= π
where γ ∈ R is a parameter that depends on the locations of the singularities of
0
R f(γ+ik∆) cosk∆t and γ = A , simplify Eq. (7.10) to give
(7.10)
(7.11)
eA/2
2t 2t t
A+2kπi
k=1
(e) Write a function (.m file) that implements Eq. 7.11. Your function should take an anonymous function, a vector of t values, and A as an input, and output a vector of f’s at t.
(f) Use your function from (e) to evaluate the inverse transforms of ˆ s−1 1
i. f(s) = (s−1)2+1 − s ˆ21
k=1
ii. f(s) = s − s+1
Compare your function with MATLAB’s inbuilt ilaplace() function.
4