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

60 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
Matlab Program 3.3 A matlab program containing the Gauss-Seidel method to solve the equations in Exercise 3.5
clear all
maxiter=1000;
A=[9 4 1; 1 6 0;
1 -2 -6];
C=[-17; 4;
14];
%%%%%%%%%%%%%%%%%%%%%%
%Initial Guess
%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(size(C));
[rows,columns]=size(A);
for k=1:maxiter
for i=1:rows
sum=0.0;
for j=1:columns
if j~=i
sum=sum-A(i,j)*x(j);
end
end
x(i)=(sum+C(i))/A(i,i);
end
residual=A*x-C
if max(abs(residual)) < 1.0e-8 fprintf(’Solution obtained !’); x break; end; end; 3.2. ITERATIVE METHODS 61 [P] = [S][Λ][S]−1 where the columns of [S] are the eigenvectors of [P] and [Λ] is a diagonal matrix whose elements are the eigenvalues of [P]. Thus we can write the error at the k’th iteration as {ε}(k+1) = [S] [Λ]k+1 [S]−1 {ε}(0) So in order for the error to go to zero as fast as possible, we would need that the magnitude of all the eigenvalues of [P ] less than 1. Exercise 3.6: Consider the system of equations in Exercise 3.5. Define the [M] and [N] matrices. Calculate the eigenvalues of [M]−1 [N] and show that they are consistent with your results in Exercise 3.5. Exercise 3.7: The equation [A] {x} = {C} (3.37) can be solved by iterative schemes which can be written as [M]{x}(k+1) =[N]{x}(k) +{C} (3.38) {x}(k) is the approximate value of the vector {x} at the k’th iteration. [A] = [M] − [N] are matrices {C} is a constant vector. 1. Consider a system of two equations where 􏰈a b􏰉 [A]=[M]−[N]= c d What is the condition among the elements of [A] for the iterative scheme to converge if 􏰈a 0􏰉 [M]= c d (3.39) 62 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS 2. Solve 􏰈17􏰉􏰊x1􏰋 􏰊4􏰋 23x=5 (3.40) 2 with the iterative scheme and the structure of matrix [M] given by Eq. (3.39). Start with 􏰊 x1 􏰋(0) 􏰊 2 􏰋 x=0 and write down the the behaviour of your approximated solution for 4 iterations. Can you get a converged solution ? 3. Swap the rows of Eq. (3.40) and try to find the solution again. Can you get a converged solution ? Explain why. 2 3.3 Newton’s method for nonlinear equations (re- visited) Earlier in the course, we talked about finding roots of two simultaneous nonlin- ear equations. What if we want to find roots of 3 or more simultaneous nonlinear equations ? The methods that we have mentioned so far can only cope with linear equations. What can we do for nonlinear equations ? For simplicity, the theory for solving a system of 3 or more simultaneous nonlinear equations will be developed for 3 × 3 system. Extension to larger systems is straight forward. Consider the following system of 3 simultaneous nonlinear equations u(x,y,z) = 0 v(x,y,z) = 0 w(x,y,z) = 0 Do a Taylor series expansion on the left hand side and truncate at the first order terms to get 3.3. NEWTON’S METHOD FOR NONLINEAR EQUATIONS (REVISITED) 63 where u(xi+1, yi+1, zi+1) v(xi+1, yi+1, zi+1) w(xi+1, yi+1, zi+1) ≈ u(xi, yi, zi) + ux(xi, yi, zi)(xi+1 − xi) + uy(xi,yi,zi)(yi+1 −yi) + uz(xi,yi,zi)(zi+1 −zi) =0 ≈ v(xi, yi, zi) + vx(xi, yi, zi)(xi+1 − xi) + vy(xi,yi,zi)(yi+1 −yi) + vz(xi,yi,zi)(zi+1 −zi) =0 ≈ w(xi, yi, zi) + wx(xi, yi, zi)(xi+1 − xi) + wy(xi,yi,zi)(yi+1 −yi) + wz(xi,yi,zi)(zi+1 −zi) =0 ux = ∂u,uy = ∂u,uz = ∂u ∂x ∂y ∂x vx = ∂v,vy = ∂v,vz = ∂v ∂x ∂y ∂z wx = ∂w,wy = ∂w,wz = ∂w ∂x ∂y ∂z We can put the above system of equations into matrix form  ux uy uz  xi+1 − xi   u   vx vy vz  yi+1 − yi = − v (3.41) wxwywz zi+1−zi w 64 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS Equation (3.41) is a linear equation so we should be able to use any of the methods discussed earlier to solve it ! Hence, to solve the original system of equations, do the following 1. Calculate all the required derivatives. 2. Guess the (initial) values for xi, yi and zi. 3. Solve Eq. (3.41) to obtain (xi+1 − xi), (yi+1 − yi) and (zi+1 − zi). You can use Gauss elimination, LU decomposition, Gauss-Seidel etc. 4. Calculate xi+1, yi+1, zi+1 since you know xi, yi and zi from step 2 above. 5. Repeat steps 2, 3 and 4 until a converged solution is obtained. This method is a more general version of the Newton-Raphson method mentioned earlier. Example 3.2: If it is desired to solve the following set of nonlinear algebraic equations 1 2 e−xy +20z+10π−3 = 0 3 using the Newton-Raphson method, one would have to define the functions u(x, y, z), v(x, y, z), w(x, y, z) and calculate the corresponding partial derivatives. Firstly, let u(x,y,z) = 3x+cos(yz)− 12 v(x,y,z) = x2 −81(y+0.1)2 +sin(z)+1.06 w(x,y,z) = e−xy+20z+10π−3. 3 The corresponding partial derivatives are 3x + cos(yz) = x2 − 81(y + 0.1)2 + sin(z) = −1.06 3.3. NEWTON’S METHOD FOR NONLINEAR EQUATIONS (REVISITED) 65 ∂u = −y sin(yz) ∂z ∂v =cos(z) ∂z ∂w = −ye−xy ∂w = −xe−xy ∂x ∂y ∂z In order to use Eq. (3.41), we would need to guess initial values of x,y and z. Starting with (x0, y0, z0) = (0, 0, 0) and substituting these numbers into Eq. (3.41) gives  3.0 0.0 0.0  x1 − x0  −0.5000000000    0.0 −16.2 1.0  y1−y0 = −0.25     0.0 0.0 20.0  z1 − z0  −10.47197551  3.0   −0.33334  0.016841 0.0046301 0.00014934  x2 − x1  0.000039     y2 − y1 =  −0.02827      z2 − z1  0.0028 ∂u=3 ∂u =−zsin(yz) ∂x∂y ∂v =2x ∂v =−16.2−162y ∂x∂y −13.464 0.16619 0.86602 20.0  3.0 0.0040504 0.00011437  x3 − x2  0.00000061    −0.33331  0.014744 −13.805 0.166246 0.86608  y3 − y2 =  0.000359     20.0  z3 − z2  −0.0000 and so on. Note that with every iteration, the right hand side becomes smaller and smaller. Eventually, the right hand side becomes a zero vector and xn+1 = xn = −0.16667, yn+1 = yn = −0.01481, zn+1 = zn = −0.523476 and we get a converged solution when n is large. ∂w = 20 Exercise 3.8: Solve the following set of nonlinear equations 66 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS xyz−x2+y2 = 1.34 xy−z2 = 0.09 ex−ey+z = 0.41 using the multi-dimensional Newton-Raphson method. One can easily extend the example above to a system of N equations. Suppose it is desired to solve f1 (x1, x2, x3, x4, ......., xN ) = 0 f2 (x1, x2, x3, x4, ......., xN ) = 0 f3 (x1,x2,x3,x4,.......,xN) = 0 , . . fN (x1, x2, x3, x4, ......., xN ) = 0 Then Eq. (3.41) extended to this system is   ∂f1 ∂f1 ∂x1 ∂x2 ∂f2 . . x2 −x2  f2 x1 ,x2 =− ,...,xN    􏰄 (k−1) (k−1) (k−1)􏰅   (k) (k−1) f x ,x ,...,x  ···112N ∂xN x1 −x1   􏰄 􏰅  . .  (k) (k−1)   (k−1) (k−1) (k−1)  ∂f1 ..  ∂x1  . . .. . .   .  .....  ∂f ··· ··· ∂f  x(k) −x(k−1) N N N N    f 􏰄 􏰅  x(k−1),x(k−1),...,x(k−1)  ∂x1 ∂xN N12N x(k−1) is the previous guess of xi and x(k) is your next (hopefully better !) guess of ii xi. The matrix of partial derivatives on the left hand side, ∂fi/∂xj, is usually known as the Jacobian. It is evaluated at x(k−1). i Chapter 4 Least Squares Suppose from an experiment you get a set of x, y values that look something like in Fig. 4.1. The data looks like it can be well approximated by a straight line. How can you draw a line a best fit through the data ? I.e. which of lines p, q or r above best fit the data ? A method of systematically calculating the line of best fit is called Least Squares approximation. y(x) q r Figure 4.1: Fitting a straight lines through a set of data points. 67 x p 68 CHAPTER 4. LEAST SQUARES 4.1 Linear least squares approximation Let’s represent any one of the lines in Fig. 4.1 as Y = ax + b (4.1) So our task now is to find a and b such that the predicted by Eq. (4.1) above is as close to the actual value of yi at a particular value of xi. The error between the predicted y value and yi is denoted by ei = yi − 􏰕􏰔􏰓􏰖 actual value (axi + b) 􏰕 􏰔􏰓 􏰖 predicted value Because we have N data points, let’s square the error at every point and sum the all up. That way we have the total error between the line and the actual data points. S = e21 + e2 + e23 + e24 + ....... + e2N N = 􏰐 e 2i i=1 N S = 􏰐 (yi − axi − b)2 i=1 We would like to find a and b, such that S is minimized. So we differentiate S with respect to a and b and the derivatives to zero. Remember that xi and yi are constants (data points) when we do the differentiation. N ∂S =􏰎2(yi−axi−b)(−xi)=0 N ∂S = 􏰎2(yi −axi −b)(−1) = 0 Exercise 4.1: Show that Eq. (4.2) can be re-written as ∂a i=1 (4.2) ∂b i=1 4.2. POLYNOMIAL LEAST SQUARES APPROXIMATION 69 NNN a 􏰐 x 2i + b 􏰐 x i = 􏰐 x i y i i=1 i=1 i=1 NN a 􏰐 xi + bN = 􏰐 yi i=1 i=1 The above two equations are just a system of 2 × 2 linear equations. They can be solved to get a and b. NNN 􏰎2􏰎 􏰎 i i􏰊a􏰋ii xxxy i=1 i=1  = i=1 􏰎N b 􏰎N i  i  xNy i=1 i=1 Exercise 4.2: Example: Suppose that you have done an experiment and collected the following data on how the resistance of a wire varies with tempera- ture T (oC) 20.5 32.7 51.0 73.2 95.7 R (ohms) 765 826 873 942 1032 Obtain a linear least squares fit of the data. 4.2 Polynomial least squares approximation Suppose we have data as shown in Fig. 4.2. Clearly, it is not appropriate to fit a straight line through the data. For this data, it is more appropriate to fit a polynomial through it. Assume a functional relationship of the form Y =a0 +a1x+a2x2 +......+anxn 70 CHAPTER 4. LEAST SQUARES 1.005 1.004 1.003 1.002 1.001 1 0.999 0.998 0.997 0.996 0.995 Figure 4.2: Data points requiring nonlinear fit of data. The error between the actual data and and the polynomial above is defined as y 0 0.2 0.4 0.6 0.8 1 x ei = yi − Yi =yi−􏰁a0+a1xi+a2x2i +......+anxni􏰂 =yi−a0−a1xi−a2x2i −......−anxni Let’s square the error and sum it NN S=􏰐e2i =􏰐􏰁yi−a0−a1xi−a2x2i −......−anxni􏰂2 i=1 i=1 (4.3) Equation (4.3) represent the total error. We would like to minimize the. To do this, calculate the partial derivatives with respect to ai’s and set them to zero. 4.2. POLYNOMIAL LEAST SQUARES APPROXIMATION 71 N ∂S =􏰎2(yi−a0−a1xi−a2x2i −......−anxni)(−1)=0, N ∂S =􏰎2(yi−a0−a1xi−a2x2i −......−anxni)(−xi)=0, N ∂S =􏰎2(yi−a0−a1xi−a2x2i −......−anxni)(−x2i)=0, . . . N ∂S =􏰎2(yi−a0−a1xi−a2x2i −......−anxni)(−xni)=0, Equation (4.4) can be written in the following matrix form ∂a0 i=1 ∂a1 i=1 ∂a2 i=1 (4.4) ∂an i=1 􏰎􏰎2􏰎3􏰎n􏰎 Nxi xi xi···xia0yi  􏰎 􏰎2 􏰎3 􏰎4 􏰎n+1  􏰎  xxxx···xaxy iiii 1ii  i 􏰎x2 􏰎x3 􏰎x4 􏰎x5 ··· 􏰎xn+2 a  􏰎x2y  iiiii2ii  􏰎x3 􏰎x4 􏰎x5 􏰎x6 ··· 􏰎xn+3  a = 􏰎x3y i i i i i3ii  . . . . . .  .   .     􏰎xn 􏰎xn+1 􏰎xn+2 􏰎xn+3 ··· 􏰎x2n  a   􏰎xny  iiiiinii The above matrix system can then be solved using techniques such as Gauss elimi- nation or LU decomposition to give the ai’s. 72 CHAPTER 4. LEAST SQUARES Chapter 5 Interpolation Sometimes, you may need to estimate the value between data points. One of the more common method is polynomial interpolation. We fit the data with the function f(x)=a0 +a1x+a2x2 +···+anxn (5.1) and require that the function passes through all the data points. See examples shown in Fig. 5.1. f(x) f(x) f(x) xxx Figure 5.1: Polynomial interpolation through a set of data points. 73 74 CHAPTER 5. INTERPOLATION 5.1 5.1.1 Figure 5.2: Linear Newton Interpolation. Newton Interpolation Linear Newton Interpolation From similar triangles, f(x1) f(x) f(x0) b x0 x x1 c a For linear interpolation, we simply draw a straight line between the two given data points a and b (see Fig. 5.2). If we want to estimate the value at point x, linear interpolation will give us the value at c whereas the true value is d. Thus there is an error associated with the linear interpolation. Shape of actual function f(x) d Linear interpolated function f1(x) f1(x) − f(x0) = f(x1) − f(x0) (5.2) x − x0 x1 − x0 Exercise 5.1: Rearrange Eq. (5.2) to give f1(x) = f(x0) + f(x1) − f(x0) (x − x0) (5.3) x1 − x0 and 5.1. NEWTON INTERPOLATION 75 􏰆 x1 − x 􏰇 􏰆 x − x0 􏰇 f1(x)=f(x0) x −x +f(x1) x −x (5.4) 10 10 A few things worth noting about the linear interpolation • Equation (5.3) is called the linear-interpolation formula • The term (f (x1 ) − f (x0 ))/(x1 − x0 ) in Eq. (5.3) is called the finite divided difference approximation of the first derivative (see later). • In general, the smaller the interval between x0 and x1, the better the approxi- mation at point x. • The subscript 1 in f1(x) denote a first order interpolating polynomial. Exercise 5.2: Use linear interpolation to estimate the value of the function f(x) = 1 − e−x at x = 1.0. Use the interval x0 = 0 and x1 = 5.0. Repeat the calculation with x1 = 4.0, x1 = 3.0 and x1 = 2.0. Illustrate on a graph how you are approaching the exact value of f(1) = 0.6321. 5.1.2 Quadratic Newton Interpolation For quadratic interpolation, let’s assume that the polynomial is of the following form f2(x)=b0 +b1(x−x0)+b2(x−x0)(x−x1) (5.5) The interpolated function is required to have the the same values as the actual function at the data points, hence • f2(x0) = f(x0) • f2(x1) = f(x1) • f2(x2) = f(x2) 76 CHAPTER 5. INTERPOLATION Shape of actual function f(x) Quadratically interpolated function f2(x) x0 x1 x2 Figure 5.3: Quadratic Newton Interpolation. Note that the subscript 2 in f2(x) denote that this is a second order interpolating polynomial. The following exercise shows that Eq. (5.5) is just another form of the polynomial function defined in Eq. (5.1). f(x2) f(x1) f(x0) Exercise 5.3: Show that Eq. (5.5) can be re-written as f2(x) = a0 + a1x + a2x2 where a0 = b0 − b1x0 + b2x0x1 a1 = b1 − b2x0 − b2x1 a2 = b2 Hence, in order to do quadratic interpolation, all we need to do is to find all the b’s in Eq. (5.5). (1) Set x = x0 in Eq. (5.5). Therefore b0 = f(x0) (5.6) 5.1. NEWTON INTERPOLATION 77 (2) If we now let x = x1 in Eq. (5.5) and use the result from step (1) above, to get b1 = f(x1) − f(x0) (5.7) x1 − x0 (3) Since f2(x2) = f(x2), we can use Eqs. (5.6) and (5.7) together with Eq. (5.5) to obtain f(x2)−f(x1) − f(x1)−f(x0) b2 = x2−x1 Equations (5.7) and (5.8) can be simplified by introducing the following notation b1 =f[x1,x0]= f(x1)−f(x0) (5.9) x1 − x0 b2 =f[x2,x1,x0]= f[x2,x1]−f[x1,x0] (5.10) x2 − x0 The advantage of using the square brackets will be clear when we look at the general form of Newton’s interpolating polynomials in the next section. 5.1.3 General form of Newton’s Interpolating Polynomials In general, if you have n data points, you will fit a polynomial of order n−1 through all the n data points. fn(x) = b0 +b1(x−x0)+b2(x−x0)(x−x1)+···+bn(x−x0)(x−x1)···(x−xn−1) where all the b’s are defined as x2 − x0 x1−x0 (5.8) where b0 = b1 = b2 = . . bn= f(x0) f[x1,x0] f[x2,x1,x0] . f[xn,xn−1,··· ,x2,x1,x0] 78 CHAPTER 5. INTERPOLATION f[xi,xj]= f(xi)−f(xj) xi − xj is called the first divided difference, f[xi,xj,xk]= f[xi,xj]−f[xj,xk] xi − xk is called the second divided difference f[xi,xj,xk,xl]= f[xi,xj,xk]−f[xj,xk,xl] xi − xl is called the third divided difference and f[xn,xn−1,...,x1,x0]= f[xn,xn−1,...,x1]−f[xn−1,xn−2,...,x0] xn − x0 is called the n’th divided difference. Notice that the second divided difference is calculated from the first divided difference. The third divided difference is calculated from the second divided difference and so on. This is illustrated in Fig. 5.4. f (xi ) First Second Third f(x0) f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0] i xi 0 x0 1 x1 2 x2 3 x3 Figure 5.4: Calculation of Divided Difference. f(x1) f[x2,x1] f[x3,x2,x1] ] f(x2) f[x3,x2] f (x3) 5.2. LAGRANGE INTERPOLATING POLYNOMIALS 79 f(x1) f(x0) x0 x1 Figure 5.5: Two data points Exercise 5.4: Use the Newton’s interpolating polynomial to fit a function through f(x) = 1 − e−x. Assume that you have data at points x =0, 1, 3 and 8. 5.2 Lagrange Interpolating Polynomials Suppose we want to interpolate between two data points, (x0,f(x0)), (x1,f(x1)) shown in Fig. 5.5. Assume that we have two functions, L0(x) and L1(x). They are defined such that 􏰊 1 if L0(x) = 0 if 􏰊 0 if L1(x) = 1 if x = x0 x = x1 x = x0 x = x1 It should be easy to convince yourself that L1(x) and L2(x) can be written as 80 CHAPTER 5. INTERPOLATION 1 L0(x) x0 L1(x) x1 Figure 5.6: Illustrative sketch of L0(x) and L1(x). 5.2. LAGRANGE INTERPOLATING POLYNOMIALS 81 f1(x)=f(x0)L0(x)+f(x1)L1(x) f(x1) f(x0)L0(x) f(x0) Figure 5.7: Illustrative sketch of Lagrange interpolation polynomial. L0(x)= (x−x1) (x0 −x1) L1(x)= (x−x0) (x1 −x0) These two functions are sketched in Fig. 5.6. The first order Lagrange interpo- lating polynomial is obtained by a linear combination of L0(x) and L1(x). f1(x) = f(x0)L0(x) + f(x1)L1(x) (5.11) A sketch of f1(x) is shown in Fig. 5.7. Suppose we three data points, (x0, f(x0)), (x1, f(x1)), (x2, f(x2)) and we want to fit a polynomial through it (see Fig. 5.8). Assume that we have three functions, L0(x),L1(x) and L2(x). They are defined such that  1 i f x = x 0 L0(x)= 0 if x=x1 (5.12) 0 if x=x2 f(x1)L1(x) x0 x1 82 CHAPTER 5. INTERPOLATION f(x2) f(x1) f(x0) x0 x1 x2 Figure 5.8: Three points for polynomial Lagrange interpolation. 0 if x=x0 L1(x)= 1 if x=x1 0 if x=x2 0 if x=x0 L2(x)= 0 if x=x1 1 if x=x2 It should be easy to convince yourself that the following and L2(x) satisfies the conditions listed in Eqs. (5.12) - (5.13) (5.14) expressions for L0(x), L1(x) (5.14). (5.15) (5.16) (5.17) in Fig. 5.9. is obtained by a linear com- L0(x) = L1(x) = L2(x) = (x − x1)(x − x2) (x0 − x1)(x0 − x2) (x − x2)(x − x0) (x1 − x2)(x1 − x0) (x − x1)(x − x0) (x2 − x1)(x2 − x0) For illustrative purposes, Eqs. (5.15) - (5.17) are shown The second order Lagrange interpolating polynomial bination of L0(x), L1(x) and L2(x). f2(x) = f(x0)L0(x) + f(x1)L1(x) + f(x2)L2(x) (5.18) 5.2. LAGRANGE INTERPOLATING POLYNOMIALS 83 1.5 11 0.5 0 0 -0.5 0 L0(x) L1(x) L2(x) xxx 012 12345678 Figure 5.9: Graphical illustration of Eqs. (5.15) - (5.17). 84 CHAPTER 5. INTERPOLATION In general, the nth order Lagrange polynomial (i.e. the Lagrange polynomial that can be fit through n + 1 data points) can be represented concisely as where n fn(x) = 􏰐 Li(x)f(xi) i=0 Li(x) = 􏰑n x − xj j=0 xi −xj j ̸=i (5.19) (5.20) Exercise 5.5: Given the following data x0=1 f(x0)=2 x1 = 4 f(x1) = 3.386294 x2 = 6 f(x2) = 3.791760 Estimate the value of f(x = 2) using both first and second order Lagrange Interpolating polynomials. The matlab program solution for this exercise is given in Matlab Program 5.1 Note that • The Lagrange interpolating polynomial is just a different form of the Newton interpolating polynomials. • The main advantage is that it is slightly easier to program. • It is slightly slower to compute than the Newton Interpolating polynomials. Exercise 5.6: Show that both the Lagrange and Newton interpolating poly- nomials are identical by reducing the first and second order Newton polynomials to first and second order Lagrange polynomials respectively. 5.3. SPLINE INTERPOLATION 85 Matlab Program 5.1 Matlab program for Exercise 5.5. clear all; xi=[1 4 6]; fxi=[2 3.386294 3.791760]; x=2 N=length(xi); fx=0; for i=1:N Lix=1.0; for j=1:N if j ~= i Lix=Lix*(x-xi(j))/(xi(i)-xi(j)); end end fx=fx+fxi(i)*Lix; end 86 CHAPTER 5. INTERPOLATION 20 15 f(x) x) 10 (f 5 Cubic spline interpolation Lagrange interpolation Exact function Data points 00 2 4 6 8 Figure 5.10: Illustration of the advantage of spline interpolation. x x 5.3. SPLINE INTERPOLATION 87 (a) polynomial function S(x) S1(x) S2(x) (b) fi(x) is a single S0(x) xx Figure 5.11: Newton or Lagrange interpolation polynomial (a) and Spline function (b). 5.3 Spline Interpolation Newton or Lagrange interpolation can lead to erroneous results when there is an abrupt change in data (see example in Fig. 5.10) Splines are made up of piecewise polynomials connecting only two data points. This is different to Newton or Lagrange polynomials where you have only one polynomial connecting all data points. See Fig. 5.11. • The spline function,S(x), is made up of several polynomial functions, Si(x). • S(x)=Si(x)forxi x(i)
j=j-1;
break end
end
f(i)=a(j)+b(j)*(x(i)-x_data(j))+c(j)*(x(i)-x_data(j))*(x(i)-x_data(j));
end
plot(x_data,f_data,’ko’,x,f);
xlabel(’x’);ylabel(’f(x)’);

96
CHAPTER 5.
INTERPOLATION
10
8
6
4
2
0
−2
−4
0 1 2 3 4 5 6 7 8 9 10
x
Figure 5.15: The output from Matlab program 5.2. quadratic spline, ◦ data points.
• Note that S(x) is continuous, the derivative of S(x), S′(x), is continuous and the second derivative of S(x), S′′(x), is also continuous at x = xi.
• If there are n + 1 data points, there are only n number of defined Si(x).
We need to find all the ai’s, bi’s, ci’s and di’s (i = 0..n−1) in order to completely determine S(x). For reasons that will become apparent later, we need to ”introduce”
f(x)

5.3. SPLINE INTERPOLATION
97
S2(x)
S0(x)
S1(x)
f(x3) f(x1)
f(x2) f(x0)
Figure 5.16: Cubic spline.
an and cn in order to completely determine S(x). Hence, there are n + 1 number of ai’s and ci’s (i = 0..n) but only n number of bi and di (i = 0..n−1) . The unknowns that need to be solve for in order to completely determine cubic splines are
a0, a1, …., an
b0, b1, …., bn−1
c0, c1, …., cn
d0, d1, …., dn−1
Note that there are n + 1 ai’s and ci’s but only n number of bi’s and di’s. Thus the total number of unknowns are (n+1)+n+(n+1)+n = 4n+2. In order to find all the unknowns, we need 4n + 2 equations. To find these equations,
1. S(x) have the values of f at x = xi, hence,
Si(xi) = ai = f(xi) (5.40) where i = 0..n. Equation (5.40) will give you n + 1 equations.
x0 x1 x2 x3

98 2.
CHAPTER 5. INTERPOLATION
the function values at adjacent polynomials must be equal at the interior knots.
Si(xi+1) = Si+1(xi+1)
ai + bi(xi+1 − xi) + ci(xi+1 − xi)2 + di(xi+1 − xi)2 = ai+1
For the sake of conciseness, let hi = xi+1 − xi, so
ai + bihi + cih2i + dih3i = ai+1 (5.41) where i = 0..n − 1. This will give us another n number of equations.
We also want the derivative of Si(x) to be continuous at the interior nodes ′′
which leads to
bi + 2cihi + 3dih2i = bi+1 (5.43) where i = 0..n−2. Note that for most cubic splines, bn does not exist. Equation
3.
4.
(5.43) will give us another n − 1 equations.
the double derivative of Si(x) to be continuous at the interior nodes.
which gives
ci + 3dihi = ci+1
where i = 0..n − 1. This will give us another n equations.
Si(xi+1) = Si+1(xi+1) (5.42)
′′ ′′
(5.44)
(5.45)
Si (xi+1) = Si+1(xi+1)
The four conditions above will give us (n+1)+n+(n−1)+n = 4n number of equations (Eqs. (5.40), (5.41), (5.43), (5.45)). We are still in need of two more equations because there is a total of 4n + 2 number of unknowns. If you know more information about the function you are trying to approximate (e.g. its derivative at the end points etc.) you can use this information to generate the two additional

5.3. SPLINE INTERPOLATION 99
information. In the absence of any information, it is common to require that the double derivative of the cubic spline is zero at the two end points. This will lead to the following two equations
c0 = 0 (5.46) cn = 0 (5.47)
These two additional equations will give us a total of 4n + 2 equations. So it is now possible to solve for all the ai’s, bi’s, ci’s and di’s.
In order to solve for all the unknowns, do the following. Re-arrange (5.45) to get
di = ci+1 − ci 3hi
Put Eq. (5.48) into Eq. (5.43) to get
bi+1 = bi + (ci + ci+1) hi Put Eq. (5.48) into Eq. (5.41) to get
bi = 1 (ai+1 −ai)− hi (2ci +ci+1) hi 3
(5.48)
(5.49)
(5.50)
Exercise 5.9: Substitute Eq. (5.50) into Eq. (5.49) and show that you will get the following relationship between ai’s and ci’s
hj−1cj−1 +2(hj +hj−1)cj +hjcj+1 = 3 (aj+1 −aj)+ 3 (aj−1 −aj) (5.51) hj hj−1
where j = 1..n − 1.
Equation (5.51) is a tridiagonal system and can be solve if we assume that c0 = cn = 0. The linear system of equation equation that we need to solve looks like the following
where
[A]{X} = {C} (5.52)

100
CHAPTER 5. INTERPOLATION
1000···0  h0 2(h0 +h1) h1 0 ··· 0 
 .. 0 h1 2(h1+h2) h2 . 0 
[A]= 0 0 … … … .  00 ···h2(h+h)h
n−2 n−2 n−1 n−1 000···01
 c 0  
   c1   . . .
{X} = .
   . . .
  c n 
0 
 3 (a − a ) + 3 (a − a )   2 1 0 1  h1 h0 
 3(a3−a2)+3(a1−a2)  h2 . h1
{C} =
Hence, to construct cubic splines, do the following
1. Make use of Eq. (5.40) and set ai = f(xi).
2. Solve Eq. (5.51) to obtain all other values of ci’s. Note that by solving Eq. (5.51) we are assuming that c0 = 0 and cn = 0. This is equivalent to saying that the second derivative at x0 and xn are zero !
3. Use Eq. (5.50) to obtain the bi’s.
4. Use Eq. (5.48) to obtain the di’s.
After steps 1-4, we can evaluate the function values at the point you are interested from the following formula Si(x) = ai + bi(x − xi) + ci(x − xi)2 + di(x − xi)3. The only thing tricky is to find the interval that corresponds to your point.
 .  
 3 (a −a )+ 3 (a −a )
 n n−1
 hn−1 0
hn−2
n−2 n−1  

5.3. SPLINE INTERPOLATION 101
Exercise 5.10: Given the set of data points below, use cubic splines to find thevalueoff atx=5.
f(xi) 3.0 2.5 4.5 1.0 7.0 2.5 9.0 0.5
xi
Exercise 5.11:
You are given the four data points shown in Fig. 5.14 and asked to fit a cubic spline through the four data points. The cubic spline can be thought as be- ing made up of three ’real’ quadratic polynomial S0(x), S1(x), S2(x) and one ’imaginary’ polynomial S3(x). The three polynomials can be written as
S0(x)=a0 +b0(x−x0)+c0(x−x0)2 +d0(x−x0)3
S1(x)=a1 +b1(x−x1)+c1(x−x1)2 +d1(x−x1)3
S2(x)=a2 +b2(x−x2)+c2(x−x2)2 +d2(x−x2)3
S3(x)=a3 +b3(x−x3)+c3(x−x3)2 +d3(x−x3)3 (a) Show that
ai =f(xi) fori=0,1,2,3
(5.53) (5.54) (5.55) (5.56)
(5.57)
(b) By requiring that adjacent Si(x) to have common values at the interior points, show that
a0 + b0h0 + c0h20 + d0h30 = a1 (5.58) a1 + b1h1 + c1h21 + d1h31 = a2 (5.59)

102 CHAPTER 5. INTERPOLATION
a2 + b2h2 + c2h2 + d2h32 = a3 (5.60) (c) By requiring that the derivative of the adjacent Si(x) to have common values
where hi = xi+1 − xi.
at the interior points, show that
b0 + 2c0h0 + 3d0h20 = b1 (5.61) b1 + 2c1h1 + 3d1h21 = b2 (5.62)
(d) By requiring that the double derivative of the adjacent Si(x) to have common values at the interior points, show that
d0 = c1 − c0 3h0
d1 = c2 − c1 3h1
b1 =b0 +(c0 +c1)h0 (5.66) b2 =b1 +(c1 +c2)h1 (5.67)
(f) Substitute Eqs. (5.63) – (5.65) into Eqs. (5.58) – (5.60) and show that
b0 = a1 − a0 − h0 (2c0 + c1) (5.68)
h0 3
b1 = a2 − a1 − h1 (2c1 + c2) (5.69) h1 3
(5.63)
(5.64)
(5.65) (e) Substitute Eqs. (5.63) – (5.65) into Eqs. (5.61) – (5.62) and show that
d2 = c3 − c2 3h2

5.3. SPLINE INTERPOLATION
103
b2 = a3 − a2 − h2 (2c2 + c3) (5.70) h2 3
(g) Use Eqs. (5.66) – (5.70) together with the assumption that c0 = 0 and c3 = 0 and show that the ci’s could be obtained by solving the following matrix equation
 1 h0
0 0 0  c0  2(h0+h1) h1 0 c1 =
 0
0 0 0 1 c3
h1 2(h1+h2) h2 c2  0
3 (a2 −a1)+ 3 (a0 −a1)  h1 h0
 h2 h1 
(5.71)
(f) Explain how you can obtain the coefficients for S0(x), S1(x), S2(x) from the results in parts (a)-(g).
3 (a3 −a2)+ 3 (a1 −a2) 0
Exercise 5.12:
Find the cubic spline that interpolates the following data points
f(xi) 1.0 2.0 3.0 4.0 4.0 9.0 7.0 2.0
xi
Exercise 5.13: A natural cubic spline, S(x), which interpolates the data given below
xi f(xi)
11 21
30

104 CHAPTER 5. INTERPOLATION
is defined to be
􏰊 S0(x)=1+P(x−1)−Q(x−1)3 ifα≤x<β S(x)= S1(x)=1+p(x−2)−34(x−2)2+q(x−2)3 ifγ≤x<λ Find the values of α, β, γ, λ, P, Q, p and q. Chapter 6 Differentiation and Integration In order to find a formula to differentiate a function, consider the Taylor series f(xi+1)=f(xi)+f′(xi)(xi+1 −xi)+ 1f′′(xi)(xi+1 −xi)2 + 1f′′′(ξ1)(xi+1 −xi)3 where 2! 3! (6.1) xi+1 = xi + ∆ f′(x)= df dx and ξ1 is somewhere in between xi and xi+1. Equation (6.1) can be re-arranged to (6.2) (6.3) obtain f ′ (xi ) = f (xi+1 ) − f (xi ) − 1 f ′′ (xi ) (xi+1 − xi ) + 1 f ′′′ (ξ1 ) (xi+1 − xi )2 (xi+1 −xi) 2! 3! Hence f′(x) can be approximated as f′(xi)= f(xi+1)−f(xi) (xi+1 − xi) Equation (6.3) is called the Forward Difference Scheme (FDS) The other terms that have been neglected in Eq. (6.2) gives an approximation of the error in approx- imating f′(xi) with Eq. (6.3). the leading term in the truncation error for the FDS approximation is 105 106 CHAPTER 6. DIFFERENTIATION AND INTEGRATION EFDS = 1f′′(ξ1)(xi+1 −xi) 2! The Taylor series expansion Eq. (6.1)could also be used to obtain an expression for f(xi−1). f(xi−1) = f(xi)+f′(xi)(xi−1 −xi) +1f′′(xi)(xi−1 −xi)2 + 1f′′′(ξ2)(xi−1 −xi)3 2! 3! = f(xi)−f′(xi)(xi −xi−1) +1f′′(xi)(xi −xi−1)2 − 1f′′′(ξ2)(xi −xi−1)3 2! 3! (6.4) Equation (6.4) can be re-arranged to obtain an another approximation for f′(xi) f′(xi)= f(xi)−f(xi−1) (6.5) (xi − xi−1) Equation (6.5) is called the Backward Difference Scheme (BDS) approximation of the first derivative. If we assume that xi−1 is very close to xi, then the leading term in the truncation error for the FDS approximation is EBDS = 1f′′(ξ2)(xi −xi−1) 2! Exercise 6.1: Subtract Eq. (6.4) from Eq. (6.1) and show that, after some algebra, that the first derivative can be written as f′(xi)= f(xi+1)−f(xi−1) −f′′(xi)(xi+1−xi)2−(xi−xi−1)2 xi+1 −xi−1 − f′′′ (ξ1) (xi+1−xi)3 − f′′′ (ξ2) 3! (xi+1 −xi−1 ) 3! 2!(xi+1 −xi−1 ) (xi−xi−1)3 (xi+1 −xi−1 ) Following Ex. (6.1), yet another formula for computing f′(xi) is f′(xi)= f(xi+1)−f(xi−1) xi+1 − xi−1 (6.6) This formula is called the Central Difference Scheme (CDS) and its leading order error is given by 6.1. APPROXIMATION OF THE 2ND DERIVATIVE 107 ECDS =f′′(xi)(xi+1 −xi)2 −(xi −xi−1)2 2!(xi+1 −xi−1) Exercise 6.2: For the case where all the xi’s are equally spaced, i.e. xi+1 − xi = ∆ = const, show that Eqs. (6.3), (6.5), and (6.6) simplifies to the following expressions for the FDS, BDS and CDS f′(xi)= f(xi+1)−f(xi) ∆ f′(xi)= f(xi)−f(xi−1) ∆ f′(xi)= f(xi+1)−f(xi−1) 2∆ Also show that the leading error term in the FDS and BDS is O(∆) and the leading error term in the CDS is O(∆2). Exercise 6.3: An alternative way of obtaining approximations to the deriva- tive of a function is to fit a Lagrange polynomial through the data points and then differentiating the resulting curve. 1. Fit a Lagrange polynomial through two data points and show that you can obtain the forward and backward differencing scheme. 2. Fit a Lagrange polynomial through three data points and show that you can obtain the central differencing scheme. 6.1 Approximation of the 2nd Derivative A formula for the 2nd derivative can be obtained by evaluating f′(x) at points halfway between xi+1 and xi and between xi and xi−1 (see Fig. 6.1). d2f df (xi+1/2) − df (xi−1/2) dx dx dx2 (xi) = 􏰁xi + 12 (xi+1 − xi)􏰂 − 􏰁xi−1 + 12 (xi − xi−1)􏰂 (6.7) 108 CHAPTER 6. DIFFERENTIATION AND INTEGRATION llmm xi-1 xi xi+1 xi-1/2 xi+1/2 Figure 6.1: Figure showing the data points used to calculate the second derivative Using Eq. (6.6) to approximate the first derivative, Eq. (6.7) can be written as d2f dx2 (xi) = f(xi+1)−f(xi) − f(xi)−f(xi−1) xi+1−xi xi−xi−1 (6.8) 12 (xi+1 − xi−1) Exercise 6.4: Show that Eq. (6.8) simplifies to d2f (xi) = f(xi+1) (xi − xi−1) − f(xi) (xi+1 − xi−1) + f(xi−1) (xi+1 − xi) dx2 12 (xi+1 − xi−1) (xi+1 − xi) (xi − xi−1) In addition, if one assumes that xi+1 − xi = xi − xi−1 = ∆, then d2f (xi) = f(xi+1) − 2f(xi) + f(xi−1) (6.9) dx2 ∆2 Exercise 6.5: Apply the forward difference scheme (FDS) twice and show that an approximation for the 2nd derivative can be written as d2f (xi) = f(xi+2) − 2f(xi+1) + f(xi) (6.10) dx2 ∆2 Similarly, apply the backward difference scheme (BDS) twice and show that an approximation for the 2nd derivative can be written as d2f (xi) = f(xi) − 2f(xi−1) + f(xi−2) (6.11) dx2 ∆2 6.2. INTEGRATION 109 f(x1) f(x0) Figure 6.2: Area under a curve f1(x) f(x1) f(x0) x0 x1 x0 x1 Area A' Figure 6.3: The trapezoidal approximation of area under a curve. 6.2 Integration Integration means finding the area under a curve between two points. Let’s call these two points x0 and x1 (see Fig. 6.2). If we only have information at these two points, then the best we can do is to fit a straight line through these two points and find the area under it. We are estimating A with A’ (see Fig. 6.3). We know that the Lagrangian interpolating polynomial for two points can be written as x0 x1 f(x1) f(x0) Area A 110 CHAPTER 6. DIFFERENTIATION AND INTEGRATION f1(x)= (x−x1)f(x0)+ (x−x0)f(x1)+1f(2)(ξ(x))(x−x0)(x−x1) (x0 −x1) (x1 −x0) 2 So to obtain the area A’ all we have to do is to integrate f1(x) from x0 to x1. 􏰒x1 x0 f(x0) 􏰒x1 f1(x)dx = (x −x ) (x−x1)dx (x−x0)dx (x − x0)(x − x1)dx 0 1 x0 f(x1) 􏰒 x1 +(x −x) 1 0 x0 1 (2) 􏰒 x1 (6.12) +2f (ξ1) x0 where ξ1 is some number in (x0,x1). Note that the Weighted Mean Value Theorem for Integrals (see page 8) have been used in evaluating the integral involving the error term. This can be done because (x−x0)(x−x1) does not change sign if x0 ≤ x ≤ x1. Exercise 6.6: Integrate Eq. (6.12) and show that 􏰒 x1 x1 − x0 ∆3 f1(x)dx = 2 (f(x0) + f(x1)) − 12 f(2)(ξ) (6.13) x0 where ∆ = x1 − x0. You might like to use the following result A+m∆ 􏰒 3􏰆m3 m2 􏰇 (x−(A+k∆))(x−(A+l∆))dx=∆ 3 − 2 (k+l)+klm A The Trapezoidal rule is derived from Eq. (6.13) and it approximates the integral of f(x) as 􏰒 x1 ∆ f(x)dx ≈ 2 (f(x0) + f(x1)) (6.14) x0 In general, at any two points, xi and xi+1, the trapezoidal rule can be written as 􏰒 xi+1 xi ∆ f(x)dx ≈ 2 (f(xi) + f(xi+1)) (6.15) 6.2. INTEGRATION 111 The error of this approximation is ∆3 12f(2)(ξ) abab Figure 6.4: Integration over larger interval. So if the size of the interval, h, is large, the error would be large as well. What happens if we want to minimize the error as we integrate over a large interval ? Consider the situation shown in Fig. 6.4. Clearly, in this situation, one cannot expect to get an accurate answer if one were to use the trapezoidal rule to integrate over the interval a ≤ x ≤ b. In order to get a more accurate answer, it is common to divide the domain into smaller intervals, usually of equal length, and apply the Trapezoidal rule in each of the subintervals (see Fig. 6.5). Area estimated by the trapezoidal rule Actual area under a curve 112 CHAPTER 6. DIFFERENTIATION AND INTEGRATION ∆1 a=x0 x1 ∆1 ∆2 ∆2 ∆2 ∆2 b=x2 a=x0 x1 x2 x3 b=x4 Figure 6.5: Decomposing a large interval into smaller intervals. Intuitively, it should be clear that as you divide the large interval up into smaller in- tervals, you should get more accurate answers than simply just applying Trapezoidal rule once over the entire large interval. What is the general formula for Trapezoidal rule if we apply it to smaller subintervals ? Consider the situation shown in Fig. 6.6. If we apply the Trapezoidal rule to each of the subintervals, A0 = ∆2 (f(x0)+f(x1)) A1 = ∆2 (f(x1)+f(x2)) A2 = ∆2 (f(x2)+f(x3)) A3 = ∆2 (f(x3)+f(x4)) The trapezoidal rule estimate of the area under the curve would be just the sum of all the Ai’s 6.3. SIMPSON’S RULE 113 A = = = = In which A0+A1+A2+A3 ∆2 (f(x0) + f(x1)) + ∆2 (f(x1) + f(x2)) + ∆2 (f(x2) + f(x3)) + ∆2 (f(x3) + f(x4)) ∆2 (f(x0) + 2f(x1) + 2f(x2) + 2f(x3) + f(x4)) ∆􏰌3􏰍 2 f(x0)+2􏰐f(xi)+f(x4) i=1 general, the formula for Trapezoidal rule applied to a large interval a < x < b is divided into N subintervals of equal size h is (6.16) b 􏰌N−1 􏰍 􏰒∆􏰐 a f(x)dx≈ 2 f(a)+2 f(xi)+f(b) i=1 Exercise 6.7: Do the following integral 10 􏰒 1dx x 2 using the Trapezoidal rule with increasing number of integrals. Compare with the exact value of 1.6094379. 6.3 Simpson’s rule The Trapezoidal rule was derived by fitting the Lagrange polynomial through two points. A more accurate method of calculating integrals would be try and fit a Lagrange polynomial through three equally spaced points. Then find the area under the Lagrange polynomial (see Fig. 6.7). What we are trying to do here is to estimate the area A under the actual function with A′ which is the area under the Lagrange polynomial (see Fig. 6.8). For three data points, the Lagrange interpolating Polynomial can be written as 114 CHAPTER 6. DIFFERENTIATION AND INTEGRATION f(x) A0 A1 A2 A3 ∆∆∆∆ a=x0 x1 x2 x3 b=x4 Figure 6.6: Decomposing a large interval into smaller intervals. The total area under the curve can be quite accurately represented by A0 + A1 + A2 + A3. f2(x) = (x−x1)(x−x2) f(x0) (x0 −x1)(x0 −x2) + (x−x0)(x−x2) f(x1) (x1 −x0)(x1 −x2) + (x−x0)(x−x1) f(x2) (x2 −x0)(x2 −x1) +61(x − x0)(x − x1)(x − x2)f(3)(ξ(x))dx 6.3. SIMPSON’S RULE 115 Actual function f(x2) f(x1) f(x0) x0 where the last term is the error term associated with the Lagrange polynomial. ∆ ∆ Approximation of the actual function with Lagrange polynomial through 3 data points. x1 x2 Figure 6.7: Fitting a Lagrange polynomial through three data points. Exercise 6.8: Assume that x0, x1 and x2 are equally spaced. Integrate between x0 and x2 and show that where x2 􏰒∆ f2(x) = 3 (f(x0) + 4f(x1) + f(x2)) x0 x2 − x1 = ∆ Ignore the error term. You might like to use the following result x1 − x0 = ∆ A+m∆ 􏰒 3􏰆m3 m2 􏰇 (x−(A+k∆))(x−(A+l∆))dx=∆ 3 − 2 (k+l)+klm A (6.17) 116 CHAPTER 6. DIFFERENTIATION AND INTEGRATION f(x2) f(x1) f(x0) f(x2) f(x1) f(x0) ∆ x0 x1 x2 Actual area A x0 x1 x2 Calculated area A’ ∆ Figure 6.8: Physical illustration of Simpsons rule. ∆ ∆ In general the simpson’s rule can be written as xi+2 􏰒∆ f2(x) = 3 (f(xi) + 4f(xi+1) + f(xi+2)) xi It is more accurate than the Trapezoidal rule. If you want to integrate a function over a large interval, it is better to split the function into smaller interval and use the Simpsons rule over smaller sub-intervals. See the example in the section on Trapezoidal rule. For example, if we divide up the domain of a function into 10 intervals, i.e. if we have 11 data points, (x0, f(x0)), (x1, f(x1)), (x2, f(x2)),....,(x10, f(x10)), (6.18) 6.3. SIMPSON’S RULE 117 􏰒 x10 f (x)dx x0 ∆ ≈ 3 + ∆3 + ∆3 + ∆3 + ∆3 (f(x0) + 4f(x1) + f(x2)) (f(x2) + 4f(x3) + f(x4)) (f(x4) + 4f(x5) + f(x6)) (f(x6) + 4f(x7) + f(x8)) (f(x8) + 4f(x9) + f(x10)) = ∆3[f(x0) + 4(f(x1)+f(x3)+f(x5)+f(x7)+f(x9)) + 2(f(x2)+f(x4)+f(x6)+f(x8)) + f(x10)] In general, the simpson’s rule can be written as  ∆ N−1 N−2 3 f(x0)+4􏰐f(xi)+2 􏰐 f(xi)+f(xN) (6.19) i=1 i=1 iodd ieven where N is the number of intervals. Note that Simpson’s rule only works of N is even i.e. the total number of data points (N + 1) is odd. Exercise 6.9: A large interval, x0 < x < xN , can be divided into N equally spaced subintervals of size ∆ with N + 1 data points (see Figure 6.9). Assuming that N is even, use Eq. (6.17) on two subintervals of size ∆ and show that the integral of f(x) over the large interval can be approximated as 􏰒 x ∆ 4∆ N − 1 N􏰐 f(x)dx ≈ 3f(x0)+ 3 f(xn) n odd 2∆ N − 2 ∆ + 􏰐 f(xn)+ f(xN) (6.20) x0 3 neven 3 118 CHAPTER 6. DIFFERENTIATION AND INTEGRATION f(xN) f(xN-1) f(x2) f(x1) f(x0) f(x) x0 x1 x2 xN-1 xN ∆ ∆ ∆ Figure 6.9: Equally spaced subintervals of size ∆ which are inside the large interval, x0 ≤x≤xN. Chapter 7 Boundary Value Problem In many engineering problems, you will have to solve ordinary differential equations that look something like where d2y + p(x)dy + q(x)y = r(x) dx2 dx (7.1) a≤x≤b We are also usually given the conditions at the boundaries y(a) = α y(b) = β Exercise 7.1: Using CDS approximation for all the derivatives, show that Eq. (7.1) can be approximated as 􏰄−pi + 2 􏰅yi−1 xi+1 −xi−1 (xi+1 −xi−1 )(xi −xi−1 ) + 􏰄− 2 (xi+1 −xi )(xi −xi−1 ) + q 􏰅 y i i (7.2) +􏰄 pi + xi+1 −xi−1 2 (xi+1 −xi−1 )(xi+1 −xi ) yi = y(xi) pi = p(xi) qi = q(xi) ri = r(xi) 􏰅yi+1 =ri where 119 120 CHAPTER 7. BOUNDARY VALUE PROBLEM From the boundary conditions, we know that y0 = α (7.3) and yN =β Eq. (7.2) together with Eq. (7.3) represent a linear system that can be solved. They can be represented by where [A]{X} = {C} 1000···0 χ1ε1η1 0···0  .. 0χ2ε2η2 .0 [A] =  0 0 ... ... ... .   0 0 · · · χ ε η  N−2 N−2 N−1 0000···1  y 0      y1   . . . {X} = .    . . .   y N   α    1  r  r  {C} = .2  .   r   N−1   β  􏰆pi 2􏰇 χi = −x −x +(x −x )(x−x ) εi = −(xi+1 −xi)(xi −xi−1) +qi 􏰆pi 2􏰇 121 i+1 i−1 i+1 i−1 i i−1 􏰆2􏰇 ηi = x −x +(x −x )(x −x) i+1 i−1 i+1 i−1 i+1 i Example 7.1: Suppose you are asked to solve d2y + 3dy + 2y = 0 (7.4) dx2 dx in the domain 0 ≤ x ≤ 1. You are also given the boundary condition y(0) = 1.0 and y(1) = e−2. From the previous chapter, you know that the derivative and double derivative can be approximated as dy(xi)= yi+1 −yi−1 (7.5) dx 2∆ d2y (xi) = yi+1 − 2yi + yi−1 (7.6) dx2 ∆2 If we now substitute Eqs. (7.5) and Eqs. (7.6) into Eq. (7.4), you will end up with 􏰆13􏰇􏰆2􏰇􏰆13􏰇 ∆2−2∆ yi−1+ −∆2+2 yi+ ∆2+2∆ yi+1=0 (7.7) Equation (7.7) is the discretised form of Eq. (7.4). Suppose we want to implement Eq. (7.7) on a grid as shown in Fig. (7.1). This grid is equally spaced with the spacing between grid points, x∆ = 0.25. The corresponding values of y on the grid shown in Fig. (7.1) are y0, y1, y2, y3, y4 respectively. The yi’s are the unknowns that you need to solve for. Thus for this grid, you have 5 unknowns. So implementing Eq. (7.7)onx=x1,x=x2,x=x3,willgive 10y0 − 30y1 + 22y2 10y1 − 30y2 + 22y3 10y2 − 30y3 + 22y4 = 0 = 0 = 0 (7.8) (7.9) (7.10) 122 CHAPTER 7. BOUNDARY VALUE PROBLEM We also know from the given boundary conditions that y0 = 1.0 y4 = e−2 You can put Eqs. (7.8)-(7.12) in matrix form as 10000y1 0      1  10−3022 0 0y0       0 10 −30 22 0  y2 = 0 0 0 10−3022y0  3  0 0 0 0 1   y     e − 2   4 One can now solve Eq. (7.13) to get y1 0 = (7.11) (7.12) (7.13)    1      y   0.5972   y2  y   3    y   4 0.3598  0.2192  (7.14) which is the approximate solution to Eq. (7.4) at x = 0, 0.25, 0.5, 0.75, 1.00.     e − 2   ∆=0.25 x0=0 x1=0.25 x2=0.5 x3=0.75 x4=1.0 Figure 7.1: Figure showing the grid used in Example 7.1 Chapter 8 Ordinary Differential Equations In many engineering problems, you will need to solve differential equations that look something like in the domain with the initial condition 8.1 Euler’s method dx = f(t,x) dt a≤t≤b x(t = a) = α (8.1) Euler’s method is probably the simplest method used to solve Eq. (8.1). Consider Taylor’s theorem x(tn+1) = x(tn) + (tn+1 − tn)dx(tn) + (tn+1 − tn)2 d2x(ξn) (8.2) dt 2dt2 where ξn is somewhere in between tn and tn+1. If we let tn+1 − tn = h and substitute Eq. (8.1) into Eq. (8.2), we get x(tn+1 ) = x(tn ) + hf (tn , xn ) + h2 d2 x (ξn ). 2 dt2 If we assume that h is small, then we can neglect the second order term in the equation above. Thus, we get the formula for Euler’s method 123 124 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS xn+1 =xn +hf(tn,xn) (8.3) where xn is the numerical approximation of the exact solution x(tn). Euler’s method is the simplest method but it is not very accurate. Matlab program 8.1 shows how Eulers method can be used to solve dx = −2x dt with the initial condition x(0) = 1.0 for 0 ≤ t ≤ 3. The matlab program compares the numerical solution with the analytical solution, which is x(t) = exp(−2t). The output from this program for h = 0.3 is shown in Fig. 8.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 00 0.5 1 1.5 2 2.5 3 t Figure 8.1: The output from Matlab program 8.1 Exercise 8.1: Using Euler’s method, solve dx = e−t dt x 8.1. EULER’S METHOD 125 Matlab Program 8.1 A matlab program illustrating the Eulers method. clear all % Sample program to do Euler integration % Parameters for simulation N=10; h=3/N; % Initial condition x(1)=1; t(1)=0; %Loop for calculating solution with Eulers method for n=1:N dxdt=-2*x(n); x(n+1)=x(n)+dxdt*h; t(n+1)=t(n)+h; end; %Specifying exact solution t_plot=linspace(0,10); x_exact=exp(-2*t_plot); %Plotting numerical solution with analytical solution plot(t_plot,x_exact,’k-’,t,x,’ko’); xlabel(’t’); ylabel(’x’); axis([0 3 0 1]); 126 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS for0 1. Hence, in order to ensure the stability of the numerical method, |σ| < 1. |σ|2 = (1+hλR)2 +(hλI)2 < 1 This is just a circle of radius 1 centred on (-1,0). This plot is called the stability diagram and is shown in Fig. 8.5 2 λIh Figure 8.5: The stability diagram of Euler method. If λ is real and negative, then in order for the numerical method to be stable, λ Rh 8.5. STABILITY AND ERROR ANALYSIS 137 h≤ 2 (8.23) |λ| Exercise 8.3: Solve with dx +2x = 0 dt (8.24) x(t = 0) = 1 Use h = 1.2, 0.8 and 0.3 in the domain 0 ≤ t ≤ 30 Exercise 8.4: Perform a stability analysis on the Crank-Nicholson method and show that it is unconditionally stable, i.e. the method is stable for all values of h. Exercise 8.5: Perform stability analysis for the second order Runge-Kutta method on the model equation dx = λx dt and show that the region of stability can be obtained by solving the following equation λ2 h2 1+λh+ 2 −eiθ =0 (8.25) Show also that for the 4th order Runge-Kutta method, the stability region is obtained by solving λ2h2 λ3h3 λ4h4 λh+ 2 + 6 + 24 +1−eiθ=0 (8.26) The solutions to Eqs. (8.25) and (8.26) are shown in Fig. 8.6. 138 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS 8.6 Systems of Ordinary Differential Equations In many engineering problems, it is essential to solve, not just one, but a set of ordinary differential equations. In general, the set of equations that you might be required to solve can be expressed as dx1 dt dx2 dt dx3 dt dx4 = f1(x1, x2, ...., xN ) = f2(x1, x2, ...., xN ) = f3(x1, x2, ...., xN ) = f4(x1, x2, ...., xN ) .=. dt .=. .=. dxN dt = fN (x1, x2, ...., xN ) If you have a linear system, then the above set of equations can be expressed in matrix-vector notation as d {x} = [A] {x} (8.27) dt where the curly brackets ({}) denotes a vector and the square brackets [] represents a matrix. {x} is the vector of dependent variables and [A] is usually a matrix of constant coefficients. Equation (8.27) can be solved using the methods described above. Applying the explicit Euler method {x}n+1 − {x}n h = [A] {x}n {x}n+1 = {x}n + h [A] {x}n At every time step, one would need to perform a matrix-vector multiplication in order to obtain the values of the vector {x} at the next time step. 8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 139 If an implicit scheme is required to solve the problem, then one would need to solve an system of linear algebraic equations at every time step. For example, applying the trapezoidal/Crank Nicholson scheme to Eq. (8.27) gives {x}n+1 − {x}n h {x}n+1 − h2 [A] {x}n+1 􏰈 h 􏰉 n+1 1 n 1 n+1 = 2[A]{x} +2[A]{x} = {x}n + h2 [A] {x}n 􏰈 h 􏰉 n I−2[A] {x} Equation (8.28) can now be solved using methods for solving a system of differ- = I+2[A] {x} ential (e.g. LU decomposition, Gauss-Seidel etc.) (8.28) Exercise 8.6: Solve the following second order ordinary differential equation with the initial conditions d2x + ω2x = 0 dt2 x(t=0) = p dx(t=0) = 0. dt (8.29) Outline the numerical algorithm for the Eulers and Trapezoidal method. 8.6.1 Multi-dimensional Runge-Kutta methods In many situations, you might be required to solve a multi-dimensional system using RK-4 and RK-2. In the examples below, it will be assumed that the system you are interested in is autonomous, i.e. the right hand side is not a function of time. The multi-dimensional form or RK-2 to solve the system of 2 equations dx1 dt dx2 dt = f1(t, x1, x2) = f2(t, x1, x2) 140 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS (i.e. extending Eq. (8.18) to a set of 2 equations) can be written as 􏰆1 1􏰇 xn+1,1 = xn,1 + 2k11 + 2k21 h 􏰆1 1􏰇 xn+1,2 = xn,2 + 2k12 + 2k22 h k11, k12, k21 and k22 can be obtained from k11 = f1(tn, xn,1, xn,2) k12 = f2(tn, xn,1, xn,2) k21 = f1 (tn +h,xn,1 +hk11,xn,2 +hk12) k22 = f2 (tn +h,xn,1 +hk11,xn,2 +hk12) Further extending the RK-2 method to solve the system of 4 equations dx1 dt dx2 dt dx3 dt dx4 dt then the solution can be written as xn+1,1 = xn+1,2 = xn+1,3 = xn+1,4 = (8.30) (8.31) = f1(t, x1, x2, x3, x4) = f2(t, x1, x2, x3, x4) = f3(t, x1, x2, x3, x4) = f4(t, x1, x2, x3, x4) xn,1 + xn,2 + xn,3 + xn,4 + 􏰆1 1􏰇 2k11 + 2k21 h 􏰆1 1􏰇 2k12 + 2k22 h 􏰆1 1􏰇 2k13 + 2k23 h 􏰆1 1􏰇 2k14 + 2k24 h (8.32) (8.33) (8.34) (8.35) 8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 141 where k11, k12, k13, k14 and k21, k22, k23 and k24 can be obtained from k21 = k22 = k23 = k24 = f1 (tn + h, xn,1 + hk11, xn,2 + hk12, xn,3 + hk13, xn,4 + hk14) f2 (tn + h, xn,1 + hk11, xn,2 + hk12, xn,3 + hk13, xn,4 + hk14) f3 (tn + h, xn,1 + hk11, xn,2 + hk12, xn,3 + hk13, xn,4 + hk14) f4 (tn + h, xn,1 + hk11, xn,2 + hk12, xn,3 + hk13, xn,4 + hk14) k11 = k12 = k13 = k14 = f1(tn, xn,1, xn,2, xn,3, xn,4) f2(tn, xn,1, xn,2, xn,3, xn,4) f3(tn, xn,1, xn,2, xn,3, xn,4) f4(tn, xn,1, xn,2, xn,3, xn,4) The RK-4 method can also be easily extended to a multi-dimensional system. Suppose you are asked to solve a system of 4 ordinary differential equations that looks like dx1 dt dx2 dt dx3 dt dx4 dt Equation (8.19) can be written in multi-dimensional form as xn+1,1 = xn+1,2 = xn+1,3 = xn+1,4 = xn,1 + xn,2 + xn,3 + xn,4 + 􏰆11 1􏰇 6k11 + 3 (k21 + k31) + 6k41 h 􏰆11 1􏰇 6k12 + 3 (k22 + k32) + 6k42 h 􏰆11 1􏰇 6k13 + 3 (k23 + k33) + 6k43 h 􏰆11 1􏰇 6k14 + 3 (k24 + k34) + 6k44 h (8.36) (8.37) (8.38) (8.39) = f1(t, x1, x2, x3, x4) = f2(t, x1, x2, x3, x4) = f3(t, x1, x2, x3, x4) = f4(t, x1, x2, x3, x4) 142 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS where the approximated solution to xi at time n + 1 is denoted as xn+1,i. The first subscript refers to the time level and the second subscript of the variable x correspond to the number of dependent variables. In order to approximate xi at various times, kij need to be calculated as follows and k21 = f1 k22 = f2 k23 = f3 k24 = f4 k31 = f1 k32 = f2 k33 = f3 k34 = f4 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k11,xn,2 + 2k12,xn,3 + 2k13,xn,4 + 2k14 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k11,xn,2 + 2k12,xn,3 + 2k13,xn,4 + 2k14 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k11,xn,2 + 2k12,xn,3 + 2k13,xn,4 + 2k14 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k11,xn,2 + 2k12,xn,3 + 2k13,xn,4 + 2k14 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k21,xn,2 + 2k22,xn,3 + 2k23,xn,4 + 2k24 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k21,xn,2 + 2k22,xn,3 + 2k23,xn,4 + 2k24 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k21,xn,2 + 2k22,xn,3 + 2k23,xn,4 + 2k24 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k21,xn,2 + 2k22,xn,3 + 2k23,xn,4 + 2k24 k11 = k12 = k13 = k14 = f1(tn, xn,1, xn,2, xn,3, xn,4) f2(tn, xn,1, xn,2, xn,3, xn,4) f3(tn, xn,1, xn,2, xn,3, xn,4) f4(tn, xn,1, xn,2, xn,3, xn,4) 8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 143 k41 = f1 k42 = f2 k43 = f3 k44 = f4 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k31,xn,2 + 2k32,xn,3 + 2k33,xn,4 + 2k34 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k31,xn,2 + 2k32,xn,3 + 2k33,xn,4 + 2k34 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k31,xn,2 + 2k32,xn,3 + 2k33,xn,4 + 2k34 􏰆hhhhh􏰇 tn + 2,xn,1 + 2k31,xn,2 + 2k32,xn,3 + 2k33,xn,4 + 2k34 Generalization of RK-4 to a system of N dependent variables should be straight forward. 144 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS Runge−Kutta 3 2 1 0 −1 −2 −3 −5 −4 −3 −2 −1 0 1 2 Figure 8.6: Figure showing the stability region of the 2nd and 4th order Runge-Kutta methods. Solutions to selected exercises 145 146 Solution to Exercise 2.7 10.0 7.5 5.0 2.5 0.0 -3 -2 -1 0 1 2 3 4 -2.5 x Solution to Exercise 2.7 147 4.0 3.5 3.0 2.5 2.0 2.0 2.5 3.0 3.5 x g_1(x) x 4.0 148 Solution to Exercise 2.7 6 4 2 0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 x -2 -4 -6 g_2(x) x Solution to Exercise 2.7 149 150 Solution to Exercise 2.7 40 30 20 10 0 2 4 6 8 10 x g_3(x) x Solution to Exercise 2.10 151 Maple worksheet to illustrate fixed point iteration scheme Define the function f ( x ) f := x/x2 K 2 x K 3 (NewtonRaphsonExample01.1) Plot f(x) to see how it looks like 10.0 7.5 5.0 2.5 0.0 -3 -2 -1 0 1 2 3 4 x -2.5 From the graph above, it is easy to see that the roots of the function f ( x ) is at x=-1 and x=3. Let's now see one can use the Newton Raphson. First we need to calculate the derivative of f(x). dfdx := x/2 x K 2 The Newton Raphson method is given by the following formula xiC1=xiK f0xi1 df 0xi1 (NewtonRaphsonExample01.2) (NewtonRaphsonExample01.3) Using this formula to iterate with the initial value x := 4 x := 3.166666667 dx (NewtonRaphsonExample01.4) 152 Solution to Exercise 2.10 x := 3.006410256 x := 3.000010240 x := 3.000000000 x := 3.000000000 x := 3.000000000 (NewtonRaphsonExample01.4) As you can clearly see, starting from x=4 will converge to the closest root at x=3. Plot graph close to the root to see if we can find out what is going on. 10.0 7.5 5.0 2.5 0.0 -2.5 2.0 2.5 3.0 3.5 x 4.0 4.5 5.0 Try to find root close to x=-2 x := K2 x := K1.166666667 x := K1.006410256 x := K1.000010240 (NewtonRaphsonExample01.5) (NewtonRaphsonExample01.6) Solution to Exercise 2.10 153 x := K1.000000000 x := K1.000000000 x := K1.000000000 (NewtonRaphsonExample01.6) As you can clearly see, starting from x=-2 will converge to the closest root at x=-1. Plot graph close to the root to see if we can find out what is going on. -3.0 -2.5 -2.0 -1.5 x -1.0 -0.5 0.0 10.0 7.5 5.0 2.5 0.0 -2.5 154 Solution to Exercise 3.3 I1 1K3L JM A := J3 1 4M (1) JM K2 0 1N Doolittle's method for LU decomposition requires that the [L] and [U] matrices be written as below for a 3x3 system Multiply the L and U matrices will give you IJ u11 temp := Jl21 u11 JKl u 31 11 u12 l21 u12 C u22 l u Cl u 31 12 32 22 u13 LM l21 u13 C u23 M (3) l u Cl u Cu MN 31 13 32 23 33 IJ 1 0 0 LM L:= Jl21 1 0 M JKl31 l32 1MN IJu11 u12 u13 LM U:= J 0 u22 u23 M (2) JK0 0 uMN 33 The product of [L] and [U] must be equal to [A]. So equate them and you will get IJ u 1 1 u 1 2 u 1 3 LM IJ 1 1 K 3 LM JluluCu luCuMJ M J2111 2112 22 J JK l u l u C l u 2113 23 M=J314M l u C l u C u M MN J JK 2 0 1 M MN 31 11 31 12 32 22 31 13 32 23 33 By comparing both sides of the above equations, you will have 9 equations with 9 unknowns. In order to solve for the 9 unknowns, you will need to follow the following sequence. Going across the first row will and compare with the [A] matrix will give you u11, u12 and u13 u11 := 1 u12 := 1 u13 := K3 (4) Going down the first column will give you l21 and l31 l21 := 3 l31 := 2 (5) Solution to Exercise 3.3 155 Going across the 2nd row will give you u22 and u23 u22 := K2 Going down the 2nd column will give you l32 Finally, comparing elements of the 3rd row and 3rd column of both sides will give ou u33 u23 := 13 (6) l32 := 1 (7) u33 := K6 (8) We have now found the [L] and [U] matrices using Doolittle's method I1 0 0L JM L= J3 1 0M JM K2 1 1N I1 1K3L JM U=J0 K2 13M (9) JM K0 0K6N Check to see if everything is OK by computing L*U-A. If LU=A, then LU-A=0. I0 0 0L JM LU K A = J0 0 0M (10) JM K0 0 0N 156 Solution to Exercise 3.4 I3L JM J2M JM J K15 M J 2 M J K5M K2N The third column of the inverse of A can be found by solving (3) IJ pK3C3qK2r LM J J J 2qCpK54 rK94K14s M M M b:= J 41617 15 5M (1) J 8 rC 8 K 2 pK 2 qC 8 sM J 1 sC17K1 pK5 qC13 rM K88228N The first column of the inverse of A can be found by solving I1 1 0 2LIaL I1L J MJMJM J2 5 2 J3 8 3 J K5 0 1 0MJbM J0M 1MJcM=J0M MJ M J M 3NKdN K0N In order to make the original equation look like the above, one would need to set p=2, q=2, r = 2 and s=1. Substitutingthesevaluesintobwillgiveyou I1L JM J1M JM J K7 M J 2 M J K1 M K2N The second column of the inverse of A can be found by solving I1 1 0 2LIaL I0L J MJMJM (2) J2 5 2 J3 8 3 J K5 0 1 0MJbM J1M 1MJcM=J0M MJ M J M 3NKdN K0N In order to make the original equation look like the above, one would need to set p=1, q=3, r = 2 and s=1. Substitutingthesevaluesintobwillgiveyou Solution to Exercise 3.4 157 I1 1 0 2LIaL I0L J MJMJM J2 5 2 0MJbM J0M J3 8 3 1MJcM=J1M J MJ M J M K5 0 1 3NKdN K0N In order to make the original equation look like the above, one would need to set p=1, q=2, r = 3 and s=1. Substitutingthesevaluesintobwillgiveyou IJ K 2 LM J K5 M J4M J 41M J 8 M J 13M K8N The fourth column of the inverse of A can be found by solving I1 1 0 2LIaL I0L J MJMJM J2 5 2 0MJbM J0M J3 8 3 1MJcM=J0M J MJ M J M K5 0 1 3NKdN K1N (4) In order to make the original equation look like the above, one would need to set p=1, q=2, r = 2 and s=2. Substitutingthesevaluesintobwillgiveyou IJ 0LM J K1 M J4M J 5M J 8M J 1M K8N Putting all the columns (Eqs. (2) - (5)) will give the inverse of A to be IJ 1 3 K2 0LM J 1 2 K5 K1M (5) J44M Ainv := JK7 K15 41 5M (6) J2 2 8 8M JK1 K5 13 1M K2288N 158 Solution to Exercise 3.4 Check to see if everything is OK by multiplying A*Ainv. This should give you the identity matrix. I1 0 0 0L JM J0 1 0 0M JM J0 0 1 0M J M K0 0 0 1N (7) Solution to Exercise 3.8 159 This worksheet outlines the steps one would take to solve the following system of equation xyz K x2 C y2 K 1.34 = 0 xy K z2 K 0.09 = 0 e x K e y C z K 0.41 = 0 Warning, the protected names norm and trace have been redefined and unprotected Warning, the name GramSchmidt has been rebound Define the functions The Jacobian, [J] is defined to be u ( x, y, z ) = x y z K x2 C y2 K 1.34 v ( x, y, z ) = x y K z2 K 0.09 w(x, y, z) =exKeyCzK0.41 IJvu vu vuLM Jvx vy vzM vyvzM Jvw vw vwM (1) [J]= Jvx Kvx vy vzN For this system [J] can be calculated to be IJyzK2x xzC2y xy LM [J] = J y x K2 zM JKex Key 1MN IJ x 0 LM IJ 1 LM The initial guess is Jy0 M = J1 M JzM JK1MN K0N Set up matrix and iterate. The system to solve at every iteration is of the form (2) Solving IJ vu 0xi,yi,zi1 J vx vu 0xi,yi,zi1 vy vu 0xi,yi,zi1 LMI L vz MJxiC1KxiM I Ju0xi,yi,zi1M vv vv MJy J MJKz L Jvx,y,z M Jvv Jvx0xi,yi,zi1 vy0xi,yi,zi1 vz0xi,yi,zi1MJiC1 KyM iM=KJ 0i i i1M Jvw0x,y,z1 vw0x,y,z1 vw0x, y,z1M iC1 Kvx i i i vy i i i vz i i i N K z MN JKw0x, y, z 1MN i i i i Jvv vv vvM 160 Solution to Exercise 3.8 gives Thus Solving gives Thus Solving K.103738540015165082 LM I J J JK 3.0425 2.4504 K2.9898 IJ x2 K x1 LM IJ Jy2Ky1M,"=", J I J J JK Jz M K2N K.7591 3.0578 JK0.950133771041813735MN 0.99275 L IJ x 3 K x 2 LM I 0 . 0 0 0 1 L I K 1 . J IJ x1 K x0 LM IJ Jy1Ky0M,"=", J J 1. JK2.7183 K2.7183 K10N 3 . 1 . L IJ x 1 K x 0 LM I 0 . 3 4 L M Jy Ky M J M 1. K2.M,J1 0M,"=",J0.09M 1.MN Jz Kz M JKK.59MN 0.0951802085692621258M Jz Kz M JKK0.0492791657229514693MN K10N IJ x1 LM IJ 0.896261459984834862 LM Jy1M,"=", J 1.09518020856926212M Jz M K1N JK0.950720834277048565MN 0.98158 L IJ x 2 K x 1 LM I 0 . 0 1 0 6 L K.7513 1.0952 0.89626 K1.9014 M, J 2 1 M, "=", J 0.01229 M 1.1003 M JyKyM J M 0.90225 K1.9003 M, J 3 2 M, "=", J K0. M M JyKyM J M 1.MN Jz Kz M JKK0.00132MN 0.00515167675034208187M Jz Kz M JKK0.000587063235234869426MN K21N IJ x2 LM IJ 0.902248052534380895 LM Jy2M,"=", J 1.10033188531960424M 2.4651 K3.0051 1.MN Jz Kz M JKK0.00013MN K32N K21N 0.00598659254954604301 LM Solution to Exercise 3.8 161 gives Thus Solving gives Thus Solving gives IJ x3 K x2 LM IJ K0.0000200216157804530296 LM Jy3Ky2M,"=", J 0.0000272899433490323116M Jz Kz M JK 0.00000136429381857177163MN K32N IJ x3 LM IJ 0.902228030918600488 LM Jy3M,"=", J 1.10035917526295335M I K.7590 J Jz M K3N 3.0580 0.90223 K3.0054 JK0.950135135335632319MN 0.99281 L IJx4Kx3LM I K0.0002L M JyKyM J M K1.9003 M, J4 3M,"=",JK0.00004M 1. MN Jz Kz M JK 0.00016MN K43N J 1.1004 JK 2.4651 I K.7590 J Jz M K4N 3.0578 JK0.950122438492582821MN 0.99271 L IJ x 5 K x 4 LM I 0 . 0 0 0 1 L IJ x4 K x3 LM Jy4Ky3M,"=", J K0.0000629366550639356396M IJ K0.00000667436618378504944 LM Jz Kz M JK K0.0000126968430495036220MN K43N IJ x4 LM IJ 0.902221356552416753 LM Jy4M,"=", J 1.10029623860788939M M JyKyM J M K54N IJx5Kx4LM IJ K0.0000104001665921580714LM Jy5Ky4M,"=", J 0.0000307533753822318448M Jz Kz M JKK0.00000194558097252621844MN K54N J 1.1003 0.90222 K1.9002 M, J 5 4 M, "=", J 0.00002 M JK 2.4651 K3.0051 1.MN Jz Kz M JKK0.00012MN 162 Solution to Exercise 3.8 Thus IJ x5 LM IJ 0.902210956385824602 LM Jy5M,"=", J 1.10032699198327166M Jz M JK0.950120492911610270MN K5N Solution to Exercise 5.12 163 Specifying the data points First define all the hi Then get all your ai Need to solve the tridiagonal system This gives Jc2M J JcM K3N IJ0LM J118M J47M J47M 3K0N Substituting into the formula for bi gives b0 := K95 Using for formula for di gives Data := [ [1, 2], [3, 4], [4, 9], [7, 2]] (1) I1 0 0 0LIJc0LM J2610MJc1M I 0L J12M J M, J M, "=", J M (4) JK22M JK0MN (5) (6) J0 1 8 3M JK0001MN JM Jc1M Jc M,"=", J2M JKcMN J M Ic0L h0 := 2 h1 := 1 h2:=3 a0 := 2 a1 := 4 a2 := 9 a3:=2 (2) (3) 141 b1 := 613 141 b2 := 535 141 JK144M 164 Solution to Exercise 5.12 d0:= 59 141 d1 := K262 141 d2 := 16 47 The cubic spline will consist of your piecewise polynomials (7) 1%x and x%3 3%x and x%4 (8) 4%x and x%7 O P 377 K 95 xC 59 (xK1)3 141 141 P 425 613 118 141 P S(x)= QK C R 2 262 3 (xK3) K (xK3) xC P K871 C 535 xK 144 (xK4)2C 16 (xK4)3 Plotting your spline and data points to see that everything is OK. P 47 141 47 141 141 141 47 47 Solution to Exercise 5.12 165 (7) 11 10 9 8 7 y6 5 4 3 2 1 1234567 x Data Cubic Spline 166 Solution to Exercise 5.13 Specifying the data points Data := [ [1, 1], [2, 1], [3, 0]] (1) Defining your piecewise polynomials "So(x)=", 1CP (xK1)CQ (xK1)3 "S1(x)=", 1Cp (xK2)C 34 (xK2)2Cq (xK2)3 (2) The cubic spline will consist of your picewise polynomials. Need to find P, Q, p and q to completely determine the cubic spline, S ( x ). O3 P 1CP (K1Cx)CQ (K1Cx) P 1%x and x%2 2%x and x%3 S(x) = Q3 (3) (4) (5) P1Cp (xK2)C R (xK2)2Cq (xK2)3 SpecifyingthatS00x11=S10x11andthatS10x21=0 gives 1CPCQ=1 74 C p C q = 0 4 Solution to Exercise 5.13 167 Specifying that S0 ' 0 x1 1 = S1 ' 0 x1 1 gives Specifying that S0 '' 0 x1 1 = S1 '' 0 x1 1 gives PC3Q=p (6) 6 Q = 32 Using Eqs. (4), (5), (6) and (7) to solve for P, Q, p and q gives abc := 4Q = 1 , P = K1 , p = 1 , q = K9 5 4424 p : = 12 q := K9 4 P := K1 4 Q : = 14 Plotting the spline and data points to make sure that everything is OK. ( 7 ) (8) ( 9 ) (10) (11) ( 1 2 ) 168 Solution to Exercise 5.13 2.0 1.8 1.6 1.4 1.2 y 1.0 0.8 0.6 0.4 0.2 0.0 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 x 2.6 2.8 3.0 Solution to Exercise 8.6 169 We are required to outline the numerical method to solve Let d2x + ω2x = 0. (8.40) dt2 x1 = x x2 = dx dt Thus Eq. 8.40 can be written as dx2 =−ω2x1. dt So the original problem can be written in terms of the new variables x1 and x2 as a system of ODE i.e. d􏰊x1􏰋 􏰈 0 1􏰉􏰊x1􏰋 = 2 (8.41) dtx2 −ω0x2 Applying Eulers method to Eq. (8.41) gives 1􏰌􏰊x1 􏰋n+1 􏰊x1 􏰋n􏰍 􏰈 0 1􏰉􏰊x1 􏰋n h x2 −=2 x2 −ω 0 x2 􏰊x1 􏰋n+1 􏰈1 0􏰉􏰊x1 􏰋n 􏰈 0 h􏰉􏰊x1 􏰋n x =01x+−ω2h0x 222 􏰊x1􏰋n+1 􏰆􏰈10􏰉􏰈0 h􏰉􏰇􏰊x1􏰋n x = 01+−ω2h0 x 22 􏰊x1􏰋n+1 􏰈 1 h􏰉􏰊x1􏰋n x =−ω2h1x (8.42) Equation (8.42) involves matrix-vector multiplication. It is given that at time t = 0, 22 􏰊 x1 􏰋0 􏰊 p 􏰋 x=0 So, using Eulers method, the estimated solution at at time t = nh can be written as 2 170 Solution to Exercise 8.6 􏰊x1􏰋n+1 􏰈 1 h􏰉n􏰊p􏰋 x =−ω2h1 0 2 If it is desired to solve Eq. (8.40) using Trapezoidal method, then 1􏰌􏰊x1 􏰋n+1 􏰊x1 􏰋n􏰍 1􏰈 0 1􏰉􏰊x1 􏰋n 1􏰈 0 1􏰉􏰊x1 􏰋n+1 􏰈 1 −h/2􏰉􏰊x1 􏰋n+1 ω2h/2 1 h x2 x2 2 −ω 0 x2 2 −ω 0 􏰈 1 h/2􏰉􏰊x1 􏰋n x2 −=2+2 = −ω2h/2 1 x 22 x 􏰊x1 􏰋n+1 1 􏰈 1 h/2􏰉􏰈 1 h/2􏰉􏰊x1 􏰋n x = 1+ω2h2/4 −ω2h/2 1 −ω2h/2 1 x 22 􏰊x1 􏰋n+1 1 􏰈1−ω2h2/2 h 􏰉􏰊x1 􏰋n x = 1+ω2h2/4 −ω2h 1−ω2h2/2 x 22 Bibliography [1] R. Burden and J. Faires. Numerical Analysis. Thomson, 1997. [2] S. Chapra and R. Canale. Numerical Methods for Engineers. McGraw-Hill, 2002. 171