Numerical Methods in Engineering (ENGR20005)
Lecture 12 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Announcement
Assignment 2
-Due 8/5/2021
-Root finding, Sys. Linear Eqn., Inerpolation, Integration
2
L11.4:
Romberg Integration: Richardson’s extrapolation
3
Richardson’s Extrapolation
The exact value of an integral can be decomposed into its numerical approximation and an error term.
I = I(Δ) + E(Δ)
(11.14)
I is the exact value of the integral
I(Δ) is the integral predicted by the integration scheme for spacing Δ E(Δ) is the error.
If we make two predictions of I, one with spacing Δ1 and also with spacing Δ2, then
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
We will assume that Δ2 < Δ1. If we look at the multiple application of the
Trapezoidal rule, we know from error analysis (see Eq. (11.8)) that
E(Δ1) ≈ Δ21 E(Δ2) Δ2
if we assume that the double derivative, f (2)(ξ) term is constant. 4
(11.15)
(11.16)
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
Substitute Eq. (11.16) into Eq. (11.15) gives
(11.15) E(Δ1) ≈ Δ21 E(Δ2) Δ2
(11.16)
I(Δ1)+E(Δ2)Δ21 ≈I(Δ2)+E(Δ2) Δ2
E(Δ2) ≈ I(Δ1) − I(Δ2) 1 − ( Δ1 )2
Substitute back into Eq. (11.15) gives
I = I(Δ2) + E(Δ2) ≈ I(Δ2) + I(Δ1) − I(Δ2) 1 − ( Δ1 )2
Δ2
Δ2
1
approximated values of the integral, I(Δ1) and I(Δ2). 5
−I(Δ1)
Trapezoidal rule. It gives you a much better estimate of the integral if you have two
I≈I(Δ2) 1+
Equation (11.17) is the Richardson extrapolation formula for the multiple application of the
( Δ1 )2 − 1 Δ2
( Δ1 )2 − 1 Δ2
(11.17)
1
Example L11.4: Use the Romberg Integration and theTrapezoidal rule to calculate integral:
I1 = ∫ 10 1x d x 2
Use Δ1 = 1.0 and Δ2 = 0.5 and compare with the exact value of 1.6094379.
I(Δ = 1.0)=1.628968253968254 I(Δ = 0.5)=1.614406323810349
From Eq. (11.17)
1+1
I ≈ I(Δ2) − I(Δ1) Δ1 =2
( Δ1 )2 − 1 Δ2
4 3
1
( Δ1 )2 − 1 Δ2
1 3
Δ2
I ≈ I(Δ = 0.5) 43 − I(Δ = 1.0) 13 = 1.609552347091047
6
L12.1: Gauss Quadrature
7
Summary of trapezoidal and Simpson’s rules.
Trapezoidal rule
∫ x 1 f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) )
x0 ΔΔ
= 2 f(x0)+ 2 f(x1)
x Simpson’s rule
f(x)dx ≈ 3 f(x0) + 4f(x1) + f(x2)
∫2 Δ( )
x0 Δ4ΔΔ
= 3 f(x0)+ 3 f(x1)+ 3 f(x2)
In general, we can write the integrals as
∫ f(x)dx≈∑wif(xi)=w0f(x0)+w1f(x1)+w2f(x2)+.... i
For the trapezoidal rule, w0,1 = Δ and for the Simpson’s rule w0,2 = Δ , w1 = 4Δ . For 233
trapezoidal and Simpson’s rules, we have also chosen (maybe because we do not know any better!) all the xis to be equally spaced separated by a constant Δ, i.e. x1 − x0 = Δ and
x2 − x0 = 2Δ etc.
What if you can freely choose xi and wi? Is there a set of optimal values of xi and wi
that can give us a good (or even exact) value of the integral?
8
The idea! ∫ab f(x)dx ≈ Δ2 (f(x0) + f(x1)) Δ = b − a
= Δ2 f(x0) + Δ2 f(x1) “Negative” error
Note that to keep things simple,
I am going to keep
w0 = w1 = Δ/2 for the upcoming
illustration. In general,
this need not be true and we are
free to choose any
values of w0 and w1
aΔb x0
x1
9
The idea!
∫ab f(x)dx ≈ Δ2 (f(x0) + f(x1)) = Δ2 f(x0) + Δ2 f(x1)
Δ = b − a “Negative” error
Δ
b
“Positive” error
“Positive” error
a
x0
x1
10
The idea!
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
“Negative” error
Δ
b
“Positive” error
“Positive” error
a
x0
x1
11
The idea!
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
“Negative” error
“Positive” error
a
Δ
“Positive” error
b
x0
x1
12
The idea!
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
“Negative” error
Δ
b
“Positive” error
“Positive” error
a
x0
x1
13
The idea!
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
“Negative” error
Δ
b
“Positive” error
“Positive” error
a
x0
x1
14
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
a
We should be able to find x0 and x1 such that the “positive” error and “negative” error
(mostly) cancels each other and result in a very good approximation of the area under
the curve (the grey shaded area).
“Negative” error
“Positive” error
“Positive” error
x0 Δ x1
b
15
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
“Negative” error
x0 Δ x1
b
“Positive” error
“Positive” error
a
We should be able to find x0 and x1 such that the “positive” error and “negative” error
(mostly) cancels each other and result in a very good approximation of the area under
the curve (the grey shaded area).
16
∫ a b f ( x ) d x ≈ Δ2 ( f ( x 0 ) + f ( x 1 ) ) = Δ2 f ( x 0 ) + Δ2 f ( x 1 )
a
We should be able to find x0 and x1 such that the “positive” error and “negative” error
(mostly) cancels each other and result in a very good approximation of the area under
curve.
“Negative” error
“Positive” error
“Positive” error
x0 Δ x1
b
17
How do we find optimum values of x0 and x1 (and w0 and w1)?
Suppose you want to calculate the integral of the function exactly by just evaluating
the function at two points and multiplying by some weight. For the purpose of this exercise
we will consider the function in the domain −1 ≤ x ≤ 1. This technique can actually be applied for calculating the integral of any domain. See later.
I = ∫1 f(x)dx = w0 f(x0) + w1f(x1) −1
We want to pick w0, w1, x0 and x1 so that the integral is exact. Note that in the previous illustration, to keep things simple, we have kept w0 = w1 = Δ/2.
In general, that need not be the case and we can find other values of w0 and w1.
We need 4 equations to solve for these 4 unknowns.
Let’s assume that f(x) can be well approximated by a polynomial
f(x) = a0 + a1x + a2x2 + a3x3 + . . . . . . ∫1 f(x)dx=∫1 a0+a1x+a2x2+a3x3+.....dx
−1 −11 1 1 1
= a0 ∫−1 1dx + a1 ∫−1 xdx + a2 ∫−1 x2dx + a3 ∫−1 x3dx + . . . . .
Since the integration is a linear operation, we can look at each of the terms individually.
18
∫1 f(x)dx=∫1 a0+a1x+a2x2+a3x3+.....dx −1 −1
=a0 +a1 +a2 +a3∫1 x3dx+..... −1
= w0 f (x0) + w1 f (x1)
Since the integration is a linear operation, we can look at each of the terms individually.
So we “pretend” that
f(x) = x3
and look at each integral independently. We are also only concerned with terms inside the integral. Hence
we can neglect all the ai in the following analysis since they are all constants that premultiplies the integral. (12.1)
(12.2) (12.3)
(12.4)
∫1 1dx −1
f(x) = 1,
f(x) = x,
f(x) = x2,
∫1 −1
xdx=0=wx +wx 0 0 1 1
∫1 xdx −1
∫1 x2dx −1
∫1 −1
1dx=2=w +w 0 1
∫1 −1
x2dx = 2 = w x2 + w x2
3
0 0 1 1
∫1 x3dx=0=w0x3+w1x3 −1 01
19
2 = w0 + w1
0 = w0x0 + w1x1
2 = w0x2 + w1x2 3 0 1
0 = w0x3 + w1x3 01
(12.1) (12.2)
(12.3) (12.4)
You can use Newton-Raphson Method to numerically solve the 4 nonlinear equations above. It is also possible to analytically solve the 4 equations above for w0, w1, x0 and x1
(see later) to get
w0 =w1 =1 x0=−x1=− 13
Hence, you can approximate the integral
∫1 f(x)dx≈f(− 1 )+f( 1 ) −1 33
(12.5)
(12.6) Equation (12.6) is called the two-point Gauss-Legendre formula and it will give you an
exact answer if f(x) is a polynomial up to third order (x3). 20
Example L12.1: Integrate the function
f(x) = 3x3 − x2 + 2
between −1 ≤ x ≤ 1 using the two-point Gauss Legendre formula and compare with Simpsons andTrapezoidal rules.
Compare also with the exact solution.
∫1 3x3−x2+2dx=10=3.333333333333 −1 3
Exact answer
Two-point Gauss-Legendre:
Simpson’s rule: Trapezoidal rule:
I≈f(−13)+f( 13)=3.33333333333
I ≈ 13 (f (−1) + 4f (0) + f (1)) = 3.33333333333333
Exact answer, requires two evaluations of the
function f(x)
Exact answer, but requires three evaluations of the
I ≈ 2 (f(−1) + f(1)) = 2 21
function f(x) Requires two
evaluations of the
function f(x), but not a good approximation.
End of Example L12.1
Example L12.2: In the previous slides, we derived the following 4 algebraic equations
2 = w0 + w1
0 = w0x0 + w1x1
23 = w 0 x 02 + w 1 x 12 0 = w 0 x 03 + w 1 x 13
that gives the optimum weights (w0 and w1) and locations (x0 and x1) to approximate the integral I=∫1 f(x)dx=w0f(x0)+w1f(x1)
−1
Use MATLAB to solve the following 4 nonlinear algebraic equations and compare with the exact solution
w0 =w1 =1 x0 = − x1 = −
1 3
23
Lecture12M02.m
w0 + w1 − 2 = 0 w0x0 + w1x1 = 0
w0 = 1 w1 = 1
x0 = − 1 x1 = 1
f(1) = x(1)+x(2)-2;
f(2) = x(1)*x(3)+x(2)*x(4); f(3)=x(1)*x(3).^2+x(2)*x(4).^2-(2./3.); f(4)=x(1)*x(3).^3+x(2)*x(4).^3;
w0x2 + w1x2 − 2 = 0 013
w0x3 + w1x3 = 0 01
xmatlab1 =
1.000000000000000e+00
1.000000000000000e+00 -5.773502700052580e-01 5.773502700052580e-01
xmatlab2 =
1.000000000000000e+00
1.000000000000000e+00 -5.773502700049091e-01 5.773502700049091e-01
f(1) = x(1)+x(2)-2;
f(2) = x(1)*x(3)+x(2)*x(4); f(3)=x(1)*x(3).^2+x(2)*x(4).^2-(2./3.); f(4)=x(1)*x(3).^3+x(2)*x(4).^3;
clear all
%asks MATLAB to print out more decimal places.
format long e
%
% It is important to have a good guess!! %
xmatlab1 = fsolve(@TwoPointGaussLegendreQuadrature,[1;1;-1;1]);
Initial guesses
options = optimoptions(@fsolve,'Display','iter','SpecifyObjectiveGradient',true);
[xmatlab2,F,exitflag,output,JAC] = fsolve(@TwoPointGaussLegendreQuadratureWithJacobian,[1;1;-1;1],options);
function [f,J] = TwoPointGaussLegendreQuadratureWithJacobian(x)
J(1,1)=1;J(1,2)=1;J(1,3)=0;J(1,4)=0; J(2,1)=x(3);J(2,2)=x(4);J(2,3)=x(1);J(2,4)=x(2); J(3,1)=x(3).^2;J(3,2)=x(4).^2;J(3,3)=2*x(1)*x(3);J(3,4)=2*x(2)*x(4); J(4,1)=x(3).^3;J(4,2)=x(4).^3;J(4,3)=3*x(1)*x(3).^2;J(4,4)=3*x(2)*x(4).^2; end
x(1) = w0 x(2) = w1 x(3) = x0 x(4) = x1
w0 =w1 =1 x0=−x1=− 13
function f = TwoPointGaussLegendreQuadrature(x)
end
24
End of Example L12.2
I have derived the two-points Gauss-Legendre Quadrature formula. It is not too difficult to extend this derivation to three-point Gauss-Legendre Quadrature formula. If you want to use three points, you would express your integral as
I = ∫1 f(x)dx = w0 f(x0) + w1f(x1) + w2 f(x2) −1
And you would need to solve for the 6 unknowns w0, w1, w2, x0, x1 and x2. The 6 equations you would need are
w0 + w1 + w2 = 2 w0x0 + w1x1 + w2x2 = 0
w0x2 + w1x2 + w2x2 = 2 0123
w0x3 + w1x3 + w2x3 = 0 0122
w0x4 + w1x4 + w2x4 = 0125
w0x5 + w1x5 + w2x5 = 0 012
Try using fsolve() to numerically solve the 6 equations above.
26
I have derived the two-points Gauss-Legendre Quadrature formula. It is not too difficult to extend this derivation to three-point Gauss-Legendre Quadrature formula. If you want to use three points, you would express your integral as
I = ∫1 f(x)dx = w0 f(x0) + w1f(x1) + w2 f(x2) −1
And you would need to solve for the equations you would need are
6 unknowns w0, w1, w2, x0, x1 and x2. The 6 w0 + w1 + w2 = 2
What if I want the data points to be at the END (-1,1)?
w0x0 + w1x1 + w2x2 = 0 w0x2 + w1x2 + w2x2 = 2
Try using fsolve() to numerically solve the 6 equations above.
0 1 2 3 w0x03 + w1x13 + w2x23 = 0
w0x4 + w1x4 + w2x4 = 2 0125
w0x5 + w1x5 + w2x5 = 0 012
27
L12.2: Legendre Polynomials
28
Legendre Polynomials
In order to get further insight about the Gauss-Legendre quadrature, I need to introduce Legendre Polynomials.
Legendre Polynomials are solution to the Legendre differential equation
(1 − x2) d2Pn(x) − 2x dPn(x) + n(n + 1)Pn(x) = 0 dx2 dx
A solution to Eq. (12.7) can be given by Rodrigues’ formula
Pn(x)= 1 dn (x2−1)n 2nn! dxn
The explicit expression for Pn(x) for small values of n is shown below P0(x) = 1
P1(x) = x
P 2 ( x ) = 12 ( 3 x 2 − 1 )
P 3 ( x ) = 12 ( 5 x 3 − 3 x )
P4(x) = 18 (35x4 − 30x2 + 3) 29
(12.7)
P1(x)
P0(x) P3(x)
P4(x)
30
P2(x)
The Legendre polynomials have the following properties
Pn(1) = 1.
Pn(x) is an n’th order polynomial.
Pn(x) has n real roots i.e. n values of x where Pn(x) = 0.
Pn(x) is an even function for even values of n and odd function for odd values of n.
P0(x) P3(x)
P4(x) P2(x)
P1(x)
31
• • • •
Let’s now look at how the Legendre Polynomials are related to Gauss Quadrature. Recall that in order to get the weights and location of the two-point Gauss-Legendre formula, we need to solve the following set of equations
w0 + w1 = 2
(12.1)
w0x0 + w1x1 = 0
(12.2)
w0x2 + w1x2 = 2
013 (12.3)
w0x3 + w1x3 = 0 01
(12.4)
Before we begin, recall that the Legendre Polynomial of order 2, P2(x) = 12 (3x2 − 1)
If we multiply Eq. (12.1) by 1/2 and then take away from Eq. (12.3) multiplied by 3/2, we get
3 (w0x2 + w1x2) − 1 (w0 + w1) = w0 1 (3x2 − 1) + w1 = 0 201220
1 (3x2 − 1)
⇒ w0P2(x0) + w1P2(x1) = 0
3 (w0x3 + w1x3) − 1 (w0x0 + w1x1) = w0x0 1 (3x2 − 1) + w1x1 = 0
(12.8) If we multiply Eq. (12.2) by 1/2 and then take away from Eq. (12.4) multiplied by 3/2, we get
2
1
1 (3x2 − 1)
2012 20
⇒ w0x0P2(x0) + w1x1P2(x1) = 0 32
(12.9)
2
1
⇒ w0P2(x0) + w1P2(x1) = 0 ⇒ w0x0P2(x0) + w1x1P2(x1) = 0
(12.8) (12.9)
You can see that if you choose x0 and x1 such that P2(x0) = P2(x1) = 0, then we will automatically satisfy Eqs. (12.8) and (12.9). Hence x0 and x1 are the roots
of the second order Legendre Polynomial P2(x).
P 2 ( x ) = 12 ( 3 x 2 − 1 ) = 0 x=± 13
x0=− 13,x1= 13
If you substitute the above values of x0 and x1 into Eq. (12.1) and (12.2), you will get
w0 =w1 =1
which is what I stated in Eq. (12.5) previously. These values of w0, w1, x0 and x1 gives
the two-point Gauss-Legendre formula
∫1 f(x)dx≈f(− 1 )+f( 1 ) −1 33
33
We can further extend this analysis for three-point Gauss-Legendre formula. For the three-point Gauss-Legendre formula, we
need to solve
w0 + w1 + w2 = 2 w0x0 + w1x1 + w2x2 = 0
w0x2 + w1x2 + w2x2 = 2 0123
w0x3 + w1x3 + w2x3 = 0 0122
w0x4 + w1x4 + w2x4 = 0125
w0x5 + w1x5 + w2x5 = 0 012
(12.10) (12.11)
(12.12) (12.13)
(12.14) (12.15)
If we multiply Eq. (12.13) by 5/2 and then take away from Eq. (12.11) multiplied by 3/2, we get 5 (w0x3 +w1x3 +w2x3)− 3 (w0x0 +w1x1 +w2x2)
20122
= w0 1 (5x3 − 3x0) + w1 1 (5x3 − 3x1) + w2 1 (5x3 − 3x2)
If we multiply Eq. (12.14) by 5/2 and then take away from Eq. (12.12) multiplied by 3/2, we get
5 (w0x4 +w1x4 +w2x4)− 3 (w0x2 +w1x2 +w2x2)
202122
=0
(12.16)
20122012 = w0x0 1 (5x3 − 3x0) + w1x1 1 (5x3 − 3x1) + w2x2 1 (5x3 − 3x2)
202122
= 0
(12.17)
34
If we multiply Eq. (12.15) by 5/2 and then take away from Eq. (12.13) multiplied by 3/2, we get 5 (w0x5 +w1x5 +w2x5)− 3 (w0x3 +w1x3 +w2x3)
20122012 = w0x2 1 (5x3 − 3x0) + w1x2 1 (5x3 − 3x1) + w2x2 1 (5x3 − 3x2)
So we would choose
020 121 222
=0
(12.18)
Recall that the third order Legendre Polynomial of order 3, P3(x) = 12 (5x3 − 3x). Substitute into Eqs. (12.16)-(12.18)
to give
w0P3(x0) + w1P3(x1) + w2P3(x2) = 0 w0x0P3(x0) + w1x1P3(x1) + w2x2P3(x2) = 0
w0x2P3(x0) + w1x2P3(x1) + w2x2P3(x2) = 0 012
P 3 ( x ) = 12 ( 5 x 3 − 3 x ) = 0
x = 0, ±
Substitute into Eqs (12.10)-(12.2) x0 = − 35 , x1 = 0,x2 = +
w 0 = w 2 = 59 35
3 5
35 to get
w 1 = 89
So the three-point Gauss-Legendre formula can be written as
1
∫−1f(x)dx≈59f(− 35)+89f(0)+59f( 35)
In general, the points of evaluation for the n-point Gauss-Legendre formula is given by the nth order Legendre polynomial. This formula is exact for polynomials up to order 2n − 1.
36
L12.3: Gauss-Legendre Quadrature for
intervals other than −1 ≤ x ≤ 1
37
Notice that we have developed the Gauss-Legendre formula for integrals between -1 and 1. This is not very useful as most integrals are for intervals more general than
−1 ≤ x ≤ 1.
We wish to use the formulae we obtained for
∫1 f(x)dx −1
∫ab g(t)dt
This can easily done be noting that the interval a ≤ t ≤ b can be transformed to
−1 ≤ x ≤ 1 using just a simple linear transformation t = px + q
We need to find p and q that relates x and t. We would like
to calculate the following integral
x=−1 when x=+1
t=a when t=b
38
x=−1 when t=a x=+1 when t=b
a=−p+q b=p+q
We can solve the above two equations to get
t = px + q Hence
Therefore
p=b−a 2
q=b+a 2
t=b−ax+b+a 22
dt = b − a dx 2
39
t = b−ax+ b+a 22
dt = b − a dx 2
Hence, in order to evaluate this integral
(12.19)
(12.20)
∫ab g(t)dt
we would use Eqs. (12.19) and (12.20) to transform the above integral into
∫bg(t)dt =∫1 g(t(x))b−adx =∫1 f(x)dx a −1 2 −1
and then we can use the two- or three-point Gauss Legendre formulae to calculate the integral of f(x) between −1 ≤ x ≤ 1
40
Example L12.3: Use Gauss quadrature to calculate the following integral.
I1 = ∫ 10 1t dt 2
We want to transform
into
∫ab g(t)dt ∫1 f(x)dx
−1
For our example, a = 2, b = 10, and g(t) = 1/t. From the previous slide,
t = b−ax+ b+a (12.19) dt = b−adx 222
Putting values of a and b into Eqs. (12.19) and (12.20) gives
(12.20)
t = 4x + 6
dt = 4dx
41
Substitute into
gives
t = 4x + 6
dt = 4dx
∫10 1dt = ∫1 2t −1
I1 = ∫ 10 1t dt 2
f (x) dx
= ∫1
∫ 1dt=∫ 4 dx=f − 1 +f 1 = 44 + 4 4 =1.565217391304348
1
4 −14x+6
4x+64dx
Two-point Gausss-Legendre
10 1
t 4x+6 ( ) ( )
2 −1 3 3−3+63+6
Three-point Gausss-Legendre
10 1
4 4x+6
d x = 5 f
9 (
5) 9
9 ( 5) 9 42
4
3
−4 5+6
+ 8 4 + 5 9−4(0)+6 9
∫ 1 dt = ∫
− 3 + 8 f (0) + 5 f 3 = 5
4 = 1.6026936026936027
2
t
−1
3
4 5+6
I1 = ∫10 1t dt 2
Error
Trapezoidal rule
2.4
0.7905620875659
Simpson’s rule
1.6888888889
0.07945097645479
Two-point Gauss Legendre
1.56521739130
0.0442205211
Three-point Gauss Legendre
1.6026936026
0.00674430974
Exact
1.6094379124341
0
Decreasing Error
43
End of Example L12.3
L12.4: Multidimensional integration
45
We want to be able to calculate the integral
I = ∫ cd ∫ ab f ( x , y ) d x d y If f (x, y) is a smooth function, then
I = ∫cd ∫ab f(x, y)dxdy = ∫ab ∫cd f(x, y)dydx
So it should not matter which integral we do first, whether the x or y integral first, you should
get the same answer.
If we perform the y integral first and let
Can do this integral using any of the integral methods.
∫cd
f(x, y)dy = g(x)
I=∫ab(∫cd f(x,y)dy)dx =∫abg(x)dx 46
Example L12.4:
The temperature distribution on a rectangular plate of length, x, 8m and width, y, 6m is given by
T(x, y) = 2xy + 2x − x2 − 2y2 + 72 Find the average temperature on the plate.
The average average temperature on the plate is given by in integral of the temperature divided by the area
A T a v e = ∫ 08 ∫ 06 T ( x , y ) d y d x Tave= 1 ∫8∫6T(x,y)dydx
6×800
The exact value of the Tave can be calculated to be Tave = 2816/48 = 58.666666667
It is useful to visualise and plot the data before you perform the integral. For illustrative purposes, we will “discretise” the temperature field at large values
of Δx and Δy. We will initially use Δx = Δy = 2. 47
The MATLAB command
will generate the following mesh
[X,Y] = meshgrid(0:2:8,0:2:6)
[X,Y]=meshgrid(0:2:8,0:2:6) X=
02468 02468 02468 02468
Y=
00000 22222 44444 66666
0
2
4
6
Δx = 2 02468x
X(1,1) Y(1,1)
X(2,1) Y(2,1)
X(3,1) Y(3,1)
X(1,2) Y(1,2)
X(2,2) Y(2,2)
(3,2) Y(3,2)
X(1,3) Y(1,3)
X(2,3) Y(2,3)
X(3,3) Y(3,3)
X(1,4) Y(1,4)
X(2,4) Y(2,4)
X(3,4) Y(3,4)
X X(4,1) X(4,2)
X(4,3) X(4,4) Y(4,3) Y(4,4)
X(1,5) Y(1,5)
X(2,5) Y(2,5)
X(3,5) Y(3,5)
X(4,5) Y(4,5)
y
NOTE that the convention in MATLAB is that increasing y in the in the direction of
increasing rows. Hence increasing y is in the downward direction. 48
Y(4,1) Y(4,2)
Δy = 2
The MATLAB command
T=2.*X.*Y+2*X-X.^2-2*Y.^2+72
will calculate the temperature T(x, y) = 2xy + 2 − x2 − 2y2 + 72 at every mesh point (the red dots).
02468x T(1,5)
0
2
4
6
y
T( T(4,1) T(4,2)
T(4,3) T(4,4)
T(2,5) T(3,5)
T(4,5)
T=2.*X.*Y+2*X-X.^2-2*Y.^2+72 T=
72 72 64 48 24
64 72 72 64 48
40 56 64 64 56
0 24 40 48 48
T(1,1)
T(2,1)
T(3,1)
T(1,2)
T(2,2)
3,2)
T(1,3)
T(2,3)
T(3,3)
T(1,4)
T(2,4)
T(3,4)
49
The MATLAB command
T=2.*X.*Y+2*X-X.^2-2*Y.^2+72
will calculate the temperature T(x, y) = 2xy + 2 − x2 − 2y2 + 72 at every mesh point (the red dots).
02468x T = 24
0
2
4
6
y
T = 0
T = 24
T = 48 T = 56
T=40 T=48 T=48
T=2.*X.*Y+2*X-X.^2-2*Y.^2+72 T=
72 72 64 48 24
64 72 72 64 48
40 56 64 64 56
0 24 40 48 48
T = 72
T = 72
T = 64
T = 48
T = 64
T = 72
T = 72
T = 64
T = 40
T = 56
T = 64
T = 64
50
The MATLAB command
will produce a contour plot of the temperature
T(x, y) = 2xy + 2 − x2 − 2y2 + 72 field.
contour(X,Y,T,'LineWidth',3)
contour(X,Y,T,'LineWidth',3) grid on
xlabel('x')
ylabel('y')
NOTE that when you plot in MATLAB, the y axis now points upwards! 51
The MATLAB command
will integrate the temperature field in the direction of the first index of T(i, j), which is i. Because i
g(x)=∫06 f(x,y)dy is in index for row (y direction), the integral is computed in the y direction. .
g=2*trapz(T)
0
2
4
02468x T = 24
T = 72
T = 72
T = 64
T = 48
T = 64
T = 72
T = 72
T = 64
T = 40
T = 56
T = 64
T = 64
g=2*trapz(T) g=
280 352 376 352 280
6
T = 48 T = 56
T=40 T=48 T=48
T = 0 ∫06 T(0,y)dy
≈ 2(72
+2(64 + 40)
+0) = 280
T = 24
y
52
∫06 T(2,y)dy ≈ 2(72
+2(72 + 56)
+24) = 352
∫06 T(4,y)dy ≈ 2(64
+2(72 + 64)
+40) = 376
∫06 T(6,y)dy ≈ 2(48
+2(64 + 64)
+48) = 352
∫06 T(8,y)dy ≈ 2(24
+2(48 + 56)
+48) = 280
Now integrate g(x) (in the x direction) and divide by the total area to give you the average temperature
Tave=dy*trapz(g)/48
will integrate the temperature field in the y direction.
Tave=1∫8g(x)dy= 1 ∫8∫6T(x,y)dydx 48 0 6×8 0 0
0
2
4
02468x T = 24
T = 72
T = 72
T = 64
T = 48
T = 64
T = 72
T = 72
T = 64
T = 40
T = 56
T = 64
T = 64
Tave=dy*trapz(g)/48 Tave =
56.6667
Final answer
Δx =Δy =2
T = 48 T = 56
T=40 T=48 T=48
2
6
T = 0 ∫06 T(0,y)dy
≈ 2(72
+2(64 + 40)
+0) =280
T = 24
86
∫0 ∫0 T(x,y)dxdy≈ 2(280+2(352+376+352)+280)=2720
∫06 T(2,y)dy ≈ 2(72
+2(72 + 56)
+24) =352
∫06 T(4,y)dy ≈ 2(64
+2(72 + 64)
+40) =376
∫06 T(6,y)dy ≈ 2(48
+2(64 + 64)
+48) =352
∫06 T(8,y)dy ≈ 2(24
+2(48 + 56)
+48) =280
y
53
End of Example L12.4