CS计算机代考程序代写 scheme matlab AI Numerical Methods in Engineering (ENGR20005)

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