Description

make sure it is:

1. uniquely named

2. m-files

3. function and script

MAE 316 Main Project

Spring 2021 Full Assignment

April 19, 2021

This is a problem about forced-convective heat transfer, as we might have in an airplane wing.

Airplane wings are generally not rectangular in shape, but this one will be rectangular to keep

things relatively simple. One end of the wing is near the engine, which exposes it to a constant

temperature, Tref . Outside, the ambient temperature is T0 , which is less than Tref . Heat is removed

from the wing via forced convection.

The flux is FourierÃ¢â‚¬â„¢s law, with a correction in the y-direction due to the changing thickness of

the wing:

q = Ã¢Ë†â€™ÃŽÂº(1 + y)Ã¢Ë†â€¡T

(1)

The source term is NewtonÃ¢â‚¬â„¢s law of cooling:

f = Ã¢Ë†â€™ÃŽÂ·(T Ã¢Ë†â€™ T0 )

(2)

where ÃŽÂ· is the convection heat transfer coefficient and T0 is the temperature of the ambient air.

Forced convection is described by a dimensionless quantity called the Nusselt number, or Nu.

The definition of the Nusselt number is a ratio of convective to conductive heat transfer in the fluid:

ÃŽÂ·y

Nu =

(3)

ÃŽÂºf

where ÃŽÂ· is the convective heat transfer coefficient, ÃŽÂºf is the heat transfer coefficient in the fluid,

and y is the length along the sheet in the direction of the fluid flow. For the particular situation of

1

forced convection over a planar surface, the Nusselt number can be approximated as a function of

the Reynolds number Re, which characterizes the flow, and the Prandtl number Pr, which depends

on the molecular properties of the fluid:

Nu = 0.332Re

1/2

1/3

Pr

= 0.332

ÃÂvy

Ã‚Âµ

1/2

Pr1/3

(4)

where ÃÂ is the fluid density, v is the velocity of the flow and Ã‚Âµ is the dynamic viscosity. Given

these two expressions for Nu along with the other flow and fluid properties, one can arrive at an

expression for the convection heat loss coefficient ÃŽÂ· in terms of y.

Derive the governing equation for the problem, including boundary conditions. Use finite differences to discretize, and write stencils for an arbitrary node in the center of the domain along with

stencils for all boundaries. Assume that the boundaries are insulated except for the left boundary,

which is fixed at Tref .

Code Implement the solution in MATLAB and produce runs for the following conditions:

Pr = 0.71

(5)

Ã‚Âµ = 1.81 Ãƒâ€” 10Ã¢Ë†â€™5

(6)

ÃÂ = 1.2041

(7)

ÃŽÂºf = 0.024

(8)

Lx = 5

(9)

Ly = 1

(10)

T0 = 0

(11)

Tref = 300

(12)

as well as for all combinations of v = {200, 400} and ÃŽÂº = {100, 200}. For the number of intervals

in the y direction use both 8 and 32 (with a regular mesh) for the run with v = 200 and ÃŽÂº = 100,

and comment on the difference you find.

2

Linear least squares is the most famous and probably the most-used curve-fitting technique. Essentially, the method enables us to find, in some sense, the Ã¢â‚¬â„¢best fitÃ¢â‚¬â„¢ of a function

to some data. That is, we find the best estimate that we can for the parameters of the

function so that the function corresponds to the data in a manner that is Ã¢â‚¬Å“optimumÃ¢â‚¬Â in the

least-squares sense (the meaning of this will be explained). The only requirement is that the

functing is linear in the parameters.

A typical example is the fitting of a polynomial to data. Polynomials are good candidates

for linear least squares because the parameters of the polynomial Ã¢â‚¬â€œ that is, the coefficients,

are linear:

f (x) = a0 + a1 x + a2 x2 + a3 x3 + Ã‚Â· Ã‚Â· Ã‚Â·

(1)

If we have some set of data Ã¢â‚¬â€œ which we might call {x, f } (the bold type indicates that

these are vectors), then fitting the above polynomial to that data entails finding some set

of coefficients (parameters) a which produces the Ã¢â‚¬Å“bestÃ¢â‚¬Â fit. The function f is linear with

respect to the parameters Ã¢â‚¬â€œ that is, there are no nonlinear combinations of parameters such

as a2 a3 .

Figure 1 shows a set of data that might be well-described Ã¢â‚¬â€œ on average Ã¢â‚¬â€œ by a straight

line. But how do we find the line that produces the Ã¢â‚¬Å“bestÃ¢â‚¬Â fit? Clearly the red line in the

figure is a better fit than the green one, but we need a metric Ã¢â‚¬â€œ a single number Ã¢â‚¬â€œ to rate

the fitness of every candidate line. We can then optimize this metric to find the line that

produces the best fit according to our criterion.

Figure 1: a dataset and two candidate fit lines

There are many ways to produce such a metric, but the most common method is that

of least squares. The optimality metric in least squares is the sum of the squares of the

1

differences between the y-values of the dataset and the y-values calcuated by the candidate

fit line at the same x-values:

N

X

Ãâ€ =

(yi Ã¢Ë†â€™ yÃŒÆ’i )2

(2)

i=1

where Ãâ€ is our metric, y is our dataset and

yÃŒÆ’ = mx + b

(3)

is our candidate line.

Finding the optimum or the minimum over all candidate parameter sets

m

ÃŽÂ¸=

b

(4)

is now our task. To accomplish this itÃ¢â‚¬â„¢s useful to introduce the problem as one of regularizing

an overdetermined system. LetÃ¢â‚¬â„¢s say that we were to naively try and set up a system of

equations to find an m and b for a line that runs through every single datapoint. This of

course is impossible, as weÃ¢â‚¬â„¢ll find out. But this system of equations would be:

y1 = mx1 + b

y2 = mx2 + b

y3 = mx3 + b

..

..

.

.

yN = mxN + b

or, in matrix-vector form,

Ã¯Â£Â¹ Ã¯Â£Â®

x1

y1

Ã¯Â£Â¯ y 2 Ã¯Â£Âº Ã¯Â£Â¯ x2

Ã¯Â£Â¯ Ã¯Â£Âº Ã¯Â£Â¯

Ã¯Â£Â¯ .. Ã¯Â£Âº = Ã¯Â£Â¯ ..

Ã¯Â£Â° . Ã¯Â£Â» Ã¯Â£Â° .

Ã¯Â£Â®

yN

xN

Ã¯Â£Â¹

1

1Ã¯Â£Âº

Ã¯Â£Âº m

.. Ã¯Â£Âº b Ã¢â€¡â€™ y = AÃŽÂ¸

.Ã¯Â£Â»

1

(5)

The system (5) is Ã¢â‚¬Å“overdeterminedÃ¢â‚¬Â because there are more equations than unknowns.

Clearly this system has no solution as written, as should be clear from the observation

that no single straight line will connect all of the datapoints in fig. 1. We require a regularization Ã¢â‚¬â€œ that is, some transformation of (5) that will reduce the size of the system

to two equations. This is what the minimization of the metric Ãâ€ provides: Ã¢â‚¬Å“least-squares

regularizationÃ¢â‚¬Â of (5).

We also recognize that we can write the metric Ãâ€ as a vector norm:

Ãâ€ = ky Ã¢Ë†â€™ yÃŒÆ’k2

2

(6)

where the notation kÃ‚Â·k indicates the Euclidean norm: the vector length. The Euclidean

norm is formed by the dot product of the vector with itself Ã¢â‚¬â€œ that is, squaring every term in

the vector and summing Ã¢â‚¬â€œ then taking the square root. The square of this norm is therefore

the summed squares of all the elements in the vector, as given in (2).

Applying that logic to the rightmost equation in (5):

Ãâ€ = (y Ã¢Ë†â€™ yÃŒÆ’)T (y Ã¢Ë†â€™ yÃŒÆ’) = (y Ã¢Ë†â€™ AÃŽÂ¸)T (y Ã¢Ë†â€™ AÃŽÂ¸) = yT y Ã¢Ë†â€™ yT AÃŽÂ¸ Ã¢Ë†â€™ ÃŽÂ¸ T AT y + ÃŽÂ¸ T AT AÃŽÂ¸

(7)

The two terms yT AÃŽÂ¸ and ÃŽÂ¸ T AT y are transposes of each other Ã¢â‚¬â€œ and thus equal, since they

are both scalars. Combining, the above equation becomes:

Ãâ€ = yT y Ã¢Ë†â€™ 2ÃŽÂ¸ T AT y + ÃŽÂ¸ T AT AÃŽÂ¸

(8)

Minimizing Ãâ€ with respect to ÃŽÂ¸ amounts to setting the partial derivatives of Ãâ€ with

respect to all of the elements of ÃŽÂ¸ equal to zero:

Ã¢Ë†â€šÃâ€

=0

Ã¢Ë†â€šÃŽÂ¸1

Ã¢Ë†â€šÃâ€

=0

Ã¢Ë†â€šÃŽÂ¸2

.. ..

. .

Ã¢Ë†â€šÃâ€

=0

Ã¢Ë†â€šÃŽÂ¸P

where P is the total number of parameters (the size of the vector ÃŽÂ¸). In matrix-vector terms

this means:

Ã¢Ë†â€™2AT y + 2AT AÃŽÂ¸ = 0 Ã¢â€¡â€™ AT AÃŽÂ¸ = AT y

(9)

Equation (9) is the formula weÃ¢â‚¬â„¢re looking for. It leads to a practical procedure for fitting

functions to data using linear least squares:

1. form the overdetermined system AÃŽÂ¸ = y

2. regularize by left-multiplying both sides by AT

3. solve the regularized system to obtain ÃŽÂ¸ = (AT A)Ã¢Ë†â€™1 AT y

3

Nonlinear problems are the most common types of problems found in engineering.

Consider the following flux q and source term g for a one-dimensional problem:

dT

dx

g = Ã¢Ë†â€™1

q = Ã¢Ë†â€™ÃŽÂºT

(1)

(2)

which is identical to the example problem used for the intro to finite differences, except that

the heat transfer coefficient now depends on the temperature. As you might expect, this

temperature-dependent behavior of the heat transfer coefficient is the rule, not the exception!

The continuity equation yields:

dT 2

d2 T

+ ÃŽÂºT 2 = 1

ÃŽÂº

(3)

dx

dx

with the corresponding second-order stencil

T Ã¢Ë†â€™ 2T + T

T Ã¢Ë†â€™ T 2

i+1

i

iÃ¢Ë†â€™1

i+1

iÃ¢Ë†â€™1

+ ÃŽÂºTi

=1

ÃŽÂº

Ã¢Ë†â€ x

Ã¢Ë†â€ x2

(4)

Equation (4) represents a system of nonlinear equations. If it were linear, it could be represented by a linear system, and we could solve the problem by implementing a matrix and

vector like we have before. But such an approach is not feasible for this system Ã¢â‚¬â€œ what

should we do?

NewtonÃ¢â‚¬â„¢s method with one unknown variable is the best way to introduce this

approach before we get to the case of solving for multiple unknows as we must in order to

find the solution to the problem represented by (4).

The method is demonstrated graphically in figure 1. The function plotted in blue is the

(nonlinear) function to be solved, f (x) = x2 , with a solution at x = 0. We start with an

initial guess, x0 = 0.8. We then draw a tangent line at the location of the initial guess.

The point at which the tangent line crosses the x-axis is the updated guess, x1 = 0.4. From

here the process starts again, with another tangent line drawn at x1 . ItÃ¢â‚¬â„¢s clear from this

illustration that eventually the process will converge at the solution.

Now that the process is intuitively clear letÃ¢â‚¬â„¢s put some math on it. We get the equation

of the tangent line by using the slope (the derivative of the function at x0 ) and a point on

the line (x0 , f (x0 )). Finding the intercept:

f (x0 ) = f 0 (x0 )x0 + b Ã¢â€¡â€™ b = f (x0 ) Ã¢Ë†â€™ f 0 (x0 )x0

(5)

yields the equation for the line:

y = f 0 (x0 )x + f (x0 ) Ã¢Ë†â€™ f 0 (x0 )x0

(6)

Setting y = 0 yields the final formula:

0 = f 0 (x0 )(x1 Ã¢Ë†â€™ x0 ) + f (x0 ) Ã¢â€¡â€™ x1 = x0 Ã¢Ë†â€™ f (x0 )/f 0 (x0 )

1

(7)

Figure 1: NewtonÃ¢â‚¬â„¢s method with one unknown

Before we move on to the multi-dimensional formula we should talk about what is really

happening here and think about a more formal way to derive the formula (7). The red line

is a linearization of the function f at x0 . The linearized version of f is easy to find the

intercept for Ã¢â‚¬â€œ we just have to solve a linear equation as reflected by (7).

How do we formally linearize a function? With a Taylor expansion:

1

f (x) = f (x0 ) + f 0 (x0 )(x Ã¢Ë†â€™ x0 ) + f 00 (x0 )(x Ã¢Ë†â€™ x0 )2 + Ã‚Â· Ã‚Â· Ã‚Â·

2

(8)

Truncating the expansion after the linear term yields:

f (x) Ã¢â€°Ë† f (x0 ) + f 0 (x0 )(x Ã¢Ë†â€™ x0 )

(9)

Setting the right-hand side of (9) equal to zero and solving for x yields the one-dimensional

formula for NewtonÃ¢â‚¬â„¢s method (7).

Multi-dimensional NewtonÃ¢â‚¬â„¢s method applies to the sets of nonlinear equations represented by (4). The first step is to arrange the equations such that they all consist of an

expression involving some (or all) of the unknowns, and setting them all equal to zero. This

means we can think of the set of equations as a vector-valued function:

f (T) = 0

2

(10)

where the function f is vector-valued in the sense that it takes a vector T as input and

returns a vector (the result of every equation) as output. The solution to the system is the

vector of temperatures T that yields the zero vector when plugged into f .

To establish an iterative method we take the same approach as in the single-variable

case: we linearize the function about an initial guess and solve the resulting linear problem,

leading to the updated guess:

f (T) Ã¢â€°Ë† f (T0 ) + f 0 (T0 )(T Ã¢Ë†â€™ T0 )

(11)

Before setting the right-hand side equal to zero and solving for T1 , we should consider the

nature of f 0 (T0 ). Clearly it is a matrix, since every term in the expression is a vector, and

the second term on the right is f 0 (T0 ) left-multiplying the vector (T Ã¢Ë†â€™ T0 ), meaning that

f 0 (T0 ) must be a matrix. It is a matrix called the Jacobian:

Ã¢Ë†â€šf1

Ã¢Ë†â€šT2

Ã¢Ë†â€šf2

Ã¢Ë†â€šT2

Ã¯Â£Â® Ã¢Ë†â€šf1

Ã¢Ë†â€šT1

Ã¢Ë†â€šf2

Ã¢Ë†â€šT1

Ã¯Â£Â¯

Ã¯Â£Â¯

f (T) = Ã¯Â£Â¯ .

Ã¯Â£Â° ..

0

Ã¢Ë†â€šfN

Ã¢Ë†â€šTN

where

Ã¯Â£Â®

Ã¯Â£Â¹

f1 (T)

Ã¯Â£Â¯ f2 (T) Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

f 0 (T) = Ã¯Â£Â¯ .. Ã¯Â£Âº ;

Ã¯Â£Â° . Ã¯Â£Â»

Ã‚Â·Ã‚Â·Ã‚Â·

Ã‚Â·Ã‚Â·Ã‚Â·

…

Ã¢Ë†â€šf1 Ã¯Â£Â¹

Ã¢Ë†â€šTN

Ã¢Ë†â€šf2 Ã¯Â£Âº

Ã¢Ë†â€šTN Ã¯Â£Âº

Ã‚Â·Ã‚Â·Ã‚Â·

Ã¢Ë†â€šfN

Ã¢Ë†â€šTN

.. Ã¯Â£Âº

. Ã¯Â£Â»

Ã¯Â£Â¹

T1

Ã¯Â£Â¯ T2 Ã¯Â£Âº

Ã¯Â£Â¯ Ã¯Â£Âº

T = Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Â° . Ã¯Â£Â»

fN (T)

(12)

Ã¯Â£Â®

(13)

TN

The form of the Jacobian comes from the multi-dimensional nature of the function f and the

nature of the Taylor expansion.

Setting the right side of (11) equal to the zero vector and re-arranging yields the multidimensional NewtonÃ¢â‚¬â„¢s method algorithm:

T1 = T0 Ã¢Ë†â€™ [f 0 (T0 )]Ã¢Ë†â€™1 f (T0 )

The way this goes in implementation is the following:

1. initialize with T0 and some tolerance Ãâ€ž

2. calculate an initial residual vector r0 = f (T0 )

3. set a residual container vector r = r0

4. WHILE krk/kr0 k > Ãâ€ž

5. compute f 0 (T0 )

6. solve f 0 (T0 )s = r

3

(14)

7. set T1 = T0 Ã¢Ë†â€™ s

8. calculate r = f (T1 )

9. shift T0 = T1

10. end WHILE

The notation kÃ‚Â·k denotes the Euclidean norm Ã¢â‚¬â€œ i.e., the length of the vector. A typical

relative tolerance Ãâ€ž is on the order of 10Ã¢Ë†â€™6 , although that can vary.

NewtonÃ¢â‚¬â„¢s method thus effectively turns the solution of a nonlinear problem into a series of

linear problems. This means that we can use all of the iterative methods for linear problems

that weÃ¢â‚¬â„¢ve learned about, so long as the Jacobian is sparse. Will it be sparse? Comparing

(4) and (12) reveals that it will be: only the partial derivatives for TiÃ¢Ë†â€™1 , Ti and Ti+1 appear

in (4), meaning that all of the other entries on the row of the Jacobian pertaining to the

particular equation represented by the stencil will be non-zero.

The other question is even more critical: is this going to work? The answer is Ã¢â‚¬ËœmaybeÃ¢â‚¬â„¢.

With NewtonÃ¢â‚¬â„¢s method there are no guarantees, except that if you start close enough to the

solution you will get there. Much of the strategy in dealing with nonlinear problems is finding

the art of getting close to the solution: in time-dependent problems, for instance, reducing

the time step gets you closer to the solution for the next time step. Time-independent

problems (those in equilibrium or steady-state) can be converted into iterative processes

employing a kind of pseudo-time, where the simulation moves stepwise from some easily

obtained solution to the more difficult target.

Algorithmically speaking, one thing that helps greatly in implementing NewtonÃ¢â‚¬â„¢s method

is ArmijoÃ¢â‚¬â„¢s rule, in which the Ã¢â‚¬Å“full Newton stepÃ¢â‚¬Â s appearing in the algorithm above is

tested to ensure that it results in a reduction of the residual over the previous step. This is

accomplished by inserting another while loop inside of the existing one above, after step 6.

The new while loop tests to make sure that the new residual is lower than the old one. If it

isnÃ¢â‚¬â„¢t, T1 is recalculated after cutting s in half. The while loop ends when the new residual

is lower or when the step has been cut to a very small length, allowing the algorithm to go

Ã¢â‚¬ËœuphillÃ¢â‚¬â„¢ a slight amount. This often can enable the algorithm to work its way out of a local

minimum.

4

Finite differences in two spatial dimensions shares a lot in common with the onedimensional version. LetÃ¢â‚¬â„¢s fix ideas right from the beginning by considering the following

flux and source expressions:

q = Ã¢Ë†â€™ÃŽÂºÃ¢Ë†â€¡T

f = 0;

(1)

(2)

Ã¢Ë†â€™Ã¢Ë†â€¡ Ã‚Â· q + f = 0

h

i

Ã¢â€¡â€™ Ã¢Ë†â€¡ Ã‚Â· Ã¢Ë†â€™ÃŽÂºÃ¢Ë†â€¡T = 0

(3)

Ã¢â€¡â€™ Ã¢Ë†â€¡2 T = 0

(5)

The continuity equation then yields:

(4)

Now imagine that this governing equation (which takes a famous form called LaplaceÃ¢â‚¬â„¢s equation) is valid on the following domain: If we want to use finite differences to solve this system

Figure 1: the domain for LaplaceÃ¢â‚¬â„¢s equation Ã¢â‚¬â€œ the origin of coordinates is at the bottom left;

the rightmost edge is a fixed temperature T1 and half the leftmost edge is a fixed temperature

T2 , whereas every other boundary is insulated

it might be easier to write it like this:

Ã¢Ë†â€š 2T

Ã¢Ë†â€š 2T

+

=0

Ã¢Ë†â€šx2

Ã¢Ë†â€šy 2

(6)

We have, of course, finite difference approximations for second derivatives, but weÃ¢â‚¬â„¢ll need to

introduce another set of indices for the second dimension. In the x-direction we can still use

the index i, whereas in the y-direction we could go for j, such that a generic node in the

center of the discretized domain and all of its neighbors looks like this: So when we discretize

(6), we get

Ti+1,j Ã¢Ë†â€™ 2Ti,j + TiÃ¢Ë†â€™1,j Ti,j+1 Ã¢Ë†â€™ 2Ti,j + Ti,jÃ¢Ë†â€™1

+

=0

(7)

Ã¢Ë†â€ x2

Ã¢Ë†â€ y 2

If we consider that we have a regular mesh, which means that Ã¢Ë†â€ x = Ã¢Ë†â€ y, we get the stencil

Ti+1,j + TiÃ¢Ë†â€™1,j Ã¢Ë†â€™ 4Ti,j + Ti,j+1 + Ti,jÃ¢Ë†â€™1 = 0

1

(8)

Figure 2: a single node on the domain, with its neighboring nodes

Implementation is a little bit different for the case of a multi-dimensional domain. The

main reason for this is the way the vector of temperatures is constructed. In the onedimensional case, we could simply create a vector out of the unknown temperatures, with

the positions of temperatures in the vector connected with the positions of nodes in the

domain:

Ã¯Â£Â®

Ã¯Â£Â¹

T1

Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯TiÃ¢Ë†â€™1 Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

T = Ã¯Â£Â¯ Ti Ã¯Â£Âº

(9)

Ã¯Â£Â¯

Ã¯Â£Âº

T

Ã¯Â£Â¯ i+1 Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â° .. Ã¯Â£Â»

Tn

Now however weÃ¢â‚¬â„¢re faced with the problem of forming a vector containing temperatures at

nodes that are arrayed in a two-dimensional mesh. How do we accomplish that?

The most reasonable method seems to be: we stack them into the vector one row at a

time:

Ã¯Â£Â®

Ã¯Â£Â¹

T1,1

Ã¯Â£Â¯ … Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯T Ã¯Â£Âº

Ã¯Â£Â¯ i,1 Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Â¯Tn,1 Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯ T1,2 Ã¯Â£Âº

T=Ã¯Â£Â¯

(10)

Ã¯Â£Âº

Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯ Ti,2 Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â°Tn,2 Ã¯Â£Â»

..

.

and so on, until we have put all rows in the domain, from the row corresponding to j = 1 to

the row corresponding to j = m, into the vector.

What does that mean for the matrix? Consider how the terms that appear in the stencil

2

(8) are arranged in the temperature vector:

Ã¯Â£Â®

..

.

Ã¯Â£Â¹

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯Ti,jÃ¢Ë†â€™1 Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯TiÃ¢Ë†â€™1,j Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯ Ti,j Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯Ti+1,j Ã¯Â£Âº

Ã¯Â£Â¯ . Ã¯Â£Âº

Ã¯Â£Â¯ .. Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯T

Ã¯Â£Âº

Ã¯Â£Â° i,j+1 Ã¯Â£Â»

..

.

(11)

How are the coefficients that appear in the stencil arranged in the matrix as a result?

Ã¯Â£Â¹

Ã¯Â£Â®

…

… … …

…

Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã¢Ë†â€™4 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â·

Ã¯Â£Âº

Ã¯Â£Â¯

Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã¢Ë†â€™4 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â·

(12)

Ã¯Â£Âº

Ã¯Â£Â¯

Ã¯Â£Âº

Ã¯Â£Â¯

Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã¢Ë†â€™4 1 Ã‚Â· Ã‚Â· Ã‚Â· 1 Ã‚Â· Ã‚Â· Ã‚Â·Ã¯Â£Â»

Ã¯Â£Â°

..

..

..

..

..

.

.

.

.

.

In other words, instead of the tri-diagonal matrix we obtained in the 1-D problem (where

every entry is zero except for the main diagonal and the two diagonals directly above &

below it), we now have two extra diagonals located n entries away Ã¢â‚¬â€œ on either side Ã¢â‚¬â€œ from

the main diagonal on every row. This is kind of an issue for Gaussian elimination, which

is O(n) complex for the tri-diagonal case but not for higher dimensions: when the nonzero

diagonal on the lower triangle is eliminated, one of the entries that was previously zero

becomes nonzero. This is the motivation for using iterative solvers that take advantage of

the matrix sparsity.

3

Finite differences in time share a lot of similarities with finite differences in space. There

are a few special considerations, however. Consider for starters the mass-spring system of

Figure 1: If we go back to the equilibrium analysis for this system, itÃ¢â‚¬â„¢s a force balance:

Figure 1: two free masses and springs Ã¢â‚¬â€œ imagine them moving

AT CAu = f

(1)

where the left-hand side of the equation is the sum of spring forces holding the masses up,

while f are the gravitational forces that pull them down. In equilibrium the forces balance,

but when they are out of balance the resultant force will create an acceleration, which can

be described by NewtonÃ¢â‚¬â„¢s second law:

M uÃŒË† = f Ã¢Ë†â€™ AT CAu

(2)

(recall that the downward direction here is positive) where the diagonal mass matrix M is:

m1 0

M=

(3)

0 m2

1

Our goal is to solve the system of equations (2) using finite differences. Since this system

is an initial-value problem, weÃ¢â‚¬â„¢ll need to first look at the nature of initial-value problems to

determine what changes we need to make in our approach.

A look at Figure 2 gives a sense of the issue we face: given u(0) Ã¢â€°Â¡ u0 , we have to come

up with a formula for u1 , using finite differences. If the problem is second-order like (2),

this implies weÃ¢â‚¬â„¢d require a second-difference approximation. For time-dependent problems,

those look like this:

un+1 Ã¢Ë†â€™ 2un + unÃ¢Ë†â€™1

(4)

uÃŒË†(tn ) Ã¢â€°Ë†

Ã¢Ë†â€ t2

and immediately it becomes clear: if n = 0 such that tn = t0 , then we require uÃ¢Ë†â€™1 , which

is a point outside of the domain. The only approximation that seems available to us is the

forward-difference approximation for a first derivative:

uÃŒâ€¡(tn ) Ã¢â€°Ë†

un+1 Ã¢Ë†â€™ un

Ã¢Ë†â€ t

(5)

which suggests that the first thing that must happen is to reduce the system to first order.

Figure 2: a grid for finite differences for initial value problems

Order reduction is possible by introducing a new time-dependent state to the system. To

illustrate, letÃ¢â‚¬â„¢s simplify the system (2) to a second-order problem with just one state. That

is, the new system will be a single mass and a single spring. Ignoring gravity, the equations

then reduce to

muÃŒË† + cu = 0

(6)

LetÃ¢â‚¬â„¢s simplify further by setting

c=m=1

(7)

uÃŒË† + u = 0

(8)

such that the equation becomes

The solution to (8) is likely to be a trig function Ã¢â‚¬â€œ both sine and cosine qualify, and the

initial conditions will determine which one we choose. The initial conditions

u(0) = 0

uÃŒâ€¡(0) = 1

(9)

(10)

u(t) = sin t

(11)

specify that the solution is

2

This analytical solution will help us assess different schemes for solving time-dependent

problems with finite differences.

To reduce the order of the problem, we introduce a new time-dependent state:

v = uÃŒâ€¡

(12)

and substitute that into (8). The new system becomes

uÃŒâ€¡ = v

vÃŒâ€¡ = Ã¢Ë†â€™u

(13)

(14)

u(0) = 0

v(0) = 1

(15)

(16)

with the initial conditions

In the Ã¢â‚¬Ëœstate spaceÃ¢â‚¬â„¢ in which we plot u(t) = sin t against v(t) = cos t, we get a perfect circle.

This stability of the exact solution Ã¢â‚¬â€œ absent any dampening or injection of energy into the

system, the mass will oscillate in a perfect sine wave forever Ã¢â‚¬â€œ is what we seek in a numerical

scheme.

Different schemes for the solution of the system of equations can be derived by choosing

different forms for the right hand side of the system (13)-(14). The derivatives, which are

on the left, are discretized with the forward difference approximation. On the right we must

therefore choose what to put in for the states: the time step n (the current time step) or

the time step n + 1 (the next time step, which weÃ¢â‚¬â„¢d like to solve for). LetÃ¢â‚¬â„¢s look at three

different options.

Forward Euler uses the current states:

un+1 Ã¢Ë†â€™ un

= vn

Ã¢Ë†â€ t

vn+1 Ã¢Ë†â€™ vn

= Ã¢Ë†â€™un

Ã¢Ë†â€ t

(17)

(18)

Re-arranging, we get:

or

un+1 = un + Ã¢Ë†â€ tvn

vn+1 = Ã¢Ë†â€™Ã¢Ë†â€ tun + vn

(19)

(20)

un+1

1

Ã¢Ë†â€ t un

=

vn+1

Ã¢Ë†â€™Ã¢Ë†â€ t 1

vn

(21)

Forward Euler is also called an explicit scheme because the formula above creates an

explicit route from the current state to the next one, without having to invert any matrices.

3

The matrix on the right-hand side weÃ¢â‚¬â„¢ll call GF . The way the scheme works is that we

multiply the current state by GF to get the next one, and then again to get the state after

that:

un+1

un

n+1 u0

2 unÃ¢Ë†â€™1

= GF

(22)

= GF

= GF

vnÃ¢Ë†â€™1

v0

vn+1

vn

ItÃ¢â‚¬â„¢s clear then that this scheme Ã¢â‚¬â€œ or any numerical scheme for a time-dependent problem Ã¢â‚¬â€œ

involves raising this matrix to an arbitrarily high power. This can tell us something about

the schemeÃ¢â‚¬â„¢s stability. Given that

GnF = SÃŽâ€ºn S Ã¢Ë†â€™1

(23)

itÃ¢â‚¬â„¢s the eigenvalues of GF that determine its stability: if any of them are greater than 1, the

scheme will add (spurious) energy into the system and it will tend to blow up; if all of the

eigenvalues are less than 1 the scheme will (spuriously) dampen the system and it will stop

oscillating Ã¢â‚¬â€œ only if all eigenvalues are exactly equal to 1 will the scheme be stable.

The eigenvalues of GF are Ã¢Ë†Å¡

easy to calculate: ÃŽÂ»F = 1 Ã‚Â± iÃ¢Ë†â€ t; these are complex numbers,

the magnitudes of which are 1 + Ã¢Ë†â€ t2 , which is greater than 1. We therefore expect the

Forward Euler scheme to be unstable: it will add numerical energy into the system and will

tend to blow up. Explicit schemes can still be used, but they must be used with caution,

with a very small timestep (the rate of the explosion can be reduced by choosing a small Ã¢Ë†â€ t)

and for short overall times.

Backward Euler is what results if we make the other choice for the right-hand side, the

state at the n + 1 timestep:

un+1 Ã¢Ë†â€™ un

= vn+1

Ã¢Ë†â€ t

vn+1 Ã¢Ë†â€™ vn

= Ã¢Ë†â€™un+1

Ã¢Ë†â€ t

(24)

(25)

Re-arranging then gives:

or

un+1 Ã¢Ë†â€™ Ã¢Ë†â€ tvn+1 = un

vn+1 + Ã¢Ë†â€ tun+1 = vn

(26)

(27)

1 Ã¢Ë†â€™Ã¢Ë†â€ t un+1

u

= n

Ã¢Ë†â€ t

1

vn+1

vn

(28)

This is also called an implicit scheme because we have to invert a matrix to find the

next state. If we call the matrix on the left above GÃ¢Ë†â€™1

B then itÃ¢â‚¬â„¢s clear that the matrix that

serves the same function as GF does for the Forward Euler is its inverse, GB . Since GÃ¢Ë†â€™1

B has

eigenvalues that are the same as those of GF , this means that all eigenvalues of GB will be

less than one, meaning that Backward Euler will be unstable with numerical dampening Ã¢â‚¬â€œ

the spring will slow down and ultimately stop as time goes to infinity.

4

Crank-Nicolson is Ã¢â‚¬â€œ literally Ã¢â‚¬â€œ an even combination of Backward and Forward Euler.

That is, half of the right-hand side of the system is states at the current time step, and half

is states at the next time step:

1

un+1 Ã¢Ë†â€™ un

= (vn + vn+1 )

Ã¢Ë†â€ t

2

vn+1 Ã¢Ë†â€™ vn

1

= Ã¢Ë†â€™ (un + un+1 )

Ã¢Ë†â€ t

2

(29)

(30)

Re-arranging:

1

1

un+1 Ã¢Ë†â€™ Ã¢Ë†â€ tvn+1 = un + Ã¢Ë†â€ tvn

2

2

1

1

Ã¢Ë†â€ tun+1 + vn+1 = Ã¢Ë†â€™ Ã¢Ë†â€ tun + vn

2

2

(31)

(32)

leading to:

1 Ã¢Ë†â€™ 12 Ã¢Ë†â€ t un+1

1

=

1

1

Ã¢Ë†â€ t

1

vn+1

Ã¢Ë†â€™ 2 Ã¢Ë†â€ t

2

1

Ã¢Ë†â€ t

2

un

1

vn

(33)

If we call the matrix on the left GCN Ã¢Ë†â€™I and the matrix on the right GCN Ã¢Ë†â€™II , then the

transfer matrix (the matrix we multiply the n time states by to get the n + 1 states) is

GÃ¢Ë†â€™1

CN Ã¢Ë†â€™I GCN Ã¢Ë†â€™II . Since GCN Ã¢Ë†â€™I and GCN Ã¢Ë†â€™II have the same eigenvalues, itÃ¢â‚¬â„¢s clear that the

eigenvalues of GÃ¢Ë†â€™1

CN Ã¢Ë†â€™I GCN Ã¢Ë†â€™II will be equal to one. Crank-Nicolson is thus the stable scheme

weÃ¢â‚¬â„¢ve been searching for.

Block matrix notation is how we write a stencil for the multi-mass problem we started

with. The first thing weÃ¢â‚¬â„¢ll do is to invoke the principle of superposition which enables us

to separate the equilibrium problem from the dynamic one: we solve the equilibrium first,

then reset the zero point for our displacements to that equilibrium point. This accounts for

the gravitational forces, so we can remove them from the dynamics given in (2). The new

dynamics then become

M uÃŒË† + AT CAu = M uÃŒË† + Ku = 0

(34)

As above we start with a reduction of the second order problem to first order by introducing

new states:

uÃŒâ€¡ = v

vÃŒâ€¡ = Ã¢Ë†â€™M Ã¢Ë†â€™1 Ku

(35)

(36)

ItÃ¢â‚¬â„¢s important to realize that if there are two masses in our system each of the vectors u and v

have two elements, meaning that each of the expressions in the system (35)-(36) as written

represent two equations (making four equations in total). Crank-Nicolson discretization

yields:

1

un+1 Ã¢Ë†â€™ un

= (vn + vn+1 )

Ã¢Ë†â€ t

2

vn+1 Ã¢Ë†â€™ vn

1 Ã¢Ë†â€™1

= Ã¢Ë†â€™ M K(un + un+1 )

Ã¢Ë†â€ t

2

5

(37)

(38)

Re-arranging:

Ã¢Ë†â€ t

Ã¢Ë†â€ t

vn+1 = un +

vn

2

2

Ã¢Ë†â€ t Ã¢Ë†â€™1

Ã¢Ë†â€ t

M Kun+1 + vn+1 = Ã¢Ë†â€™ M Ã¢Ë†â€™1 Kun + vn

2

2

un+1 Ã¢Ë†â€™

which finally yields

I

Ã¢Ë†â€™ Ã¢Ë†â€ t

I

2

Ã¢Ë†â€ t

Ã¢Ë†â€™1

M K

I

2

un+1

I

=

Ã¢Ë†â€ t

vn+1

Ã¢Ë†â€™ 2 M Ã¢Ë†â€™1 K

Ã¢Ë†â€ t

I

2

un

vn

I

(39)

(40)

(41)

Equation (41) features stacked vectors and block matrices: each written entry in the

vectors is itself a 2 Ãƒâ€” 1 vector, and each entry appearing in the matrices is a 2 Ãƒâ€” 2 matrix.

Thus the entire system is 4 Ãƒâ€” 4, although itÃ¢â‚¬â„¢s written as if it were a 2 Ãƒâ€” 2. The shorthand is

perhaps not necessary in this case, but it would be if there were many more than two springs

in the system. This is a stencil for the masses-and-springs dynamic system, which can be

used for coding an implementation of the Crank-Nicolson scheme.

6

Purchase answer to see full

attachment