MATLAB PROBLEM

1 Background

Linear algebraic equations can arise when modeling distributed systems. Linear algebraic equations can arise when modeling distributed systems.

Figure 1: A noninsulated uniform rod positioned between two walls of constant but different temperature. The finite-difference representation employs four interior nodes.

For example, Fig. 1 shows a long, thin rod positioned between two walls that are held at constant tempera- tures. Heat flows through the rod as well as between the rod and the surrounding air. For the steady-state case, a differential equation based on heat conservation can be written for such a system as

d2T +h(Ta T)=0 (1)dx2

where T is the temperature in C, x is the distance along the rod in m, h is the heat transfer coefficient between the rod and the surrounding air in m1), and Ta is the air temperature in C.

Given values for the parameters, forcing functions, and boundary conditions, calculus can be used to develop an analytical solution. For example, for h = 0.01, Ta = 20, T (0) = 40, and T (10) = 200,

T = 73.4523 exp (0.1x) 53.4523 exp (0.1x) + 20 (2)

Although it provided a solution here, calculus does not work for all such problems. In such instances, numerical methods provide a valuable alternative. In this case study, we will use finite differences to transform this differential equation into a tridiagonal system of linear algebraic equations which can be readily solved using numerical methods.

Equation (1) can be transformed into a set of linear algebraic equations by conceptualizing the rod as consisting of a series of nodes. For example, the rod in Fig. 1 is divided into six equispaced nodes. Since the rod has a length of 10, the spacing between nodes is x = 2. Calculus was necessary to solve Eq. (1) because it includes a second derivative. Finite-difference approximations provide a means to transform derivatives into algebraic form. For example, the second derivative at each node can be approximated as1

d2T = Ti+1 2Ti+Ti1 (3)

dx2

x2

1See Appendix A

1

where Ti designates the temperature at node i. This approximation can be substituted into Eq. (1) to giveTi+1 2Ti+Ti1 +h(Ta Ti)=0 (4)

x2Collecting terms and substituting the parameters gives

Ti1 + 2.04Ti Ti+1 = 0.8 (5) Thus, Eq. (1) has been transformed from a differential equation into an algebraic equation. Equation (3)

can now be applied to each of the interior nodes:

T0 + 2.04T1 T2 = 0.8 (6)T1 + 2.04T2 T3 = 0.8 (7)T2 + 2.04T3 T4 = 0.8 (8)T3 + 2.04T4 T5 = 0.8 (9)

The values of the fixed end temperatures, T0 = 40 and T5 = 200, can be substituted and moved to the right-hand side. The results are four equations with four unknowns expressed in matrix form as

2.041 0

0

1 0 2.04 1

1 2.04 0 1

0 T1 40.8

0 T2= 0.8 (10)1T3 0.8

2.04 T4 200.8

Therefore, the original differential equation has been converted into an equivalent system of linear alge- braic equations. Consequently, we can use the techniques already described to solve for the temperatures.

In addition to being a linear system, notice that Eq. (10) is also tridiagonal. The system is tridiagonal because each node depends only on its adjacent nodes. Such cases often occur when solving differential equations based on conservation laws.

2 Assignment

The following differential equation results from a force balance for a beam with uniform loading,

d2y wLx wx2 0=EIdx2 2 + 2 (11)

where x is distance along the beam in m, y is the deflection in m, L is the length in m, E is the modulus of elasticity in N m1, I is the moment of inertia in m4, and w is a uniform load in N m1.

Using Matlab

1. (a) Convert this differential equation to an equivalent system of simultaneous algebraic equations using a centered difference approximation for the second derivative.

2. (b) Develop a function to solve these equations from x = 0 to L and return the resulting distances and deflections. The first line of your function should be

3. (c) Develop a script that invokes this function and then plots the results.

4. (d) Test your script for the following parameters: L = 3m, x = 0.2m, E = 250109 Nm1, I =

   3 104 m4, w = 22 500 N m1, y(0) = 0 m, and y(3) = 0 m. 2

beam.m

function [x,y] = beam(E,I,w,y0,yL,L,dx)

A Appendix

The Taylor series expansion of a function f(x) was presented when we discussed Truncation Errors. The general expression for the expansion of f(xi+1) on the basis of the value of f(x) at xi is given as

f (xi+1 ) = f (xi ) + f (xi )h + f (xi ) h2 + f (3) (xi ) h3 + (12)2! 3!

where h xi+1 xi is called the step size. Truncating the series after the first derivative term yields

f(xi+1) = f(xi) + f(xi)h + O(h2) (13)

where the the nomenclature O(hn) represents a truncation error of the order of hn. The latter expression provides an approximation to the first derivative. Rearranging Eq. (13) we have

f(xi) = f(xi+1) f(xi) + O(h) (14)h

which is termed a forward difference because it utilizes data from i and i + 1 to estimate the derivative. A similar expression is obtained by a Taylor series expansion for xi1, i.e., a backward difference from

which yields

f (xi1 ) = f (xi ) f (xi )h + f (xi ) h2 f (3) (xi ) h3 + (15)2! 3!

f(xi) = f(xi) f(xi1) + O(h) (16)h

An approximation to the second derivative can be found by adding Eqs. (12) and (15) and rearranging the result as

f(xi) = f(xi+1) 2f(xi) + f(xi1) + O(h2) (17)h2