Source Files

Laboratory 7

Slides for this presentation are available here: 7.FiniteElements.ppt, 8.PoissonsEquation.ppt and 9.TentTestFunctions.ppt.

Interpolation in One Variable

Suppose we wish to interpolate $n$ points $(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)$. One way to find such an interpolating function is to find $n$ different functions $L_1(x), L_2(x), \ldots, L_n(x)$ where the function $L_k(x)$ passes through the point $(x_k, y_k)$ and is zero at the remaining $n - 1$ points.

This is not so hard as it sounds: suppose our we are interpolating three points $(1, 7.5), (2, 3.8), (4, 9.6)$. Consider the polynomial $\pi_1(x) = (x - 2)(x - 4)$ which is zero at both $x = 2$ and $x = 4$. If we now define $L_1(x) = 7.5 \frac{\pi_1(x)}{\pi_1(1)} = 7.5 \frac{(x - 2)(x - 4)}{(1 - 2)(1 - 4)} = 2.5(x - 2)(x - 4)$, we see that $L_1(1) = 7.5$ while $L_1(2) = L_1(4) = 0$.

Similarly, we can find define $\pi_2(x) = (x - 1)(x - 4)$ and then $L_2(x) = 3.8 \frac{\pi_2(x)}{\pi_2(2)} = 3.8 \frac{(x - 1)(x - 4)}{(2 - 1)(2 - 4)} = -1.9(x - 1)(x - 4)$ and finally $\pi_3(x) = (x - 1)(x - 2)$ and then $L_3(x) = 9.6 \frac{\pi_3(x)}{\pi_3(4)} = 9.6 \frac{(x - 1)(x - 2)}{(4 - 1)(4 - 2)} = 1.6(x - 1)(x - 2)$.

Thus, the interpolating polynomial is $p(x) = L_1(x) + L_2(x) + L_3(x) = 2.5(x - 2)(x - 4) - 1.9(x - 1)(x - 4) + 1.6(x - 1)(x - 2)$.

This technique (called Lagrange polynomial interpolation) is usually only useful if we are interpolating only two points or if the ordinate is zero for many of the points.

Interpolating two points is simpler: Consider interpolating $(2, 3.8), (4, 9.6)$. This is simply $3.8 \frac{x - 2}{1 - 2} + 9.6 \frac{x - 4}{2 - 4}$. There are a few special cases: if one of the ordinate values is 0, for example, $(2, 0), (4, 9.6)$. This is simply $0 \cdot \frac{x - 2}{1 - 2} + 9.6 \frac{x - 4}{2 - 4} = 9.6 \frac{x - 4}{2 - 4}$. If the other ordinate is 1, for example, $(2, 0), (4, 1)$, the interpolating straight line is $0 \cdot \frac{x - 2}{1 - 2} + 1 \cdot \frac{x - 4}{2 - 4} = \frac{x - 4}{2 - 4}$.

Finite Elements in One Dimensions

All of the finite-difference methods we have looked at previously include grids which are equally spaced. For simplification, we maintained equal spacing $h$ in all space dimensions; however, this is not necessary. As we have discussed previously in Laboratory 2, the domain over which we are solving the differential equation may not suit a grid and consequently, civil and aeronautical engineers came up with a more appropriate solution: the method of finite elements. For this technique, we partition the space not by a grid but by:

Intervals of arbitrary length in one dimension,
Triangular regions in two dimensions, and
Tetrahedral regions in three dimensions.

Dividing an entire region into such regions is generally called:

A partition in one dimension,
A triangular tessellation in two dimensions, and
A tetrahedral tessellation in three dimensions.

Another justification for using tessellations (even in one dimension) is that they also allow the user to select for a higher density of points (and therefore greater accuracy) at regions where we may expect the solution to change rapidly (for example, at boundaries between where two materials meet). Examples are shown in Figures 1 and 2.

Figure 1. A copper/aluminium bar between two heat sources with a partition which is concentrated at the ends and the boundary.

Figure 2. A triangular tessellation where the concentration of points is at the boundary of the two media.

This week, we will look at finite elements in one dimension; that is, a partition which divides an interval. Previously, we used equally spaced points and this yielded a very simple formula for the second derivative:

$\frac{d^2 u}{dx^2} \approx \frac{u(x + h) - 2u(x) + u(x - h)}{h^2}$

where $h$ is the distance between the points. However, if the points are not regularly spaced, the formula is slightly more complicated:

$\frac{d^2 u}{dx^2} \approx 2\frac { \left( x_k - x_{k-1} \right) u \left( x_{k+1} \right) - \left( x_{k+1} - x_{k-1} \right) u \left( x_k \right) + \left( x_{k+1} - x_k \right) u \left( x_{k-1} \right) }{ \left( x_k-x_{k-1} \right) \left( x_{k+1}-x_{k- 1} \right) \left( x_{k+1}-x_k \right) }$

(1)

With a little effort, it is quite easy to see that this formula simplifies into the previous formula if the points are equally spaced. However, such general formulas explode in higher dimensions. (You may click here to see the general formula approximating the Laplacian at $u(x_1, y_1) = z_1$ and if you were to replace $(x_2, y_2), (x_3, y_3), (x_4, y_4), (x_5, y_5)$ with $(x_1 + h, y_1), (x_1 - h, y_1), (x_1, y_1 + h), (x_1 + h, y_1 - h)$, you do get the formula $\frac{-4z_1 + z_2 + z_3 + z_4 + z_5}{h^2}$, as expected.) The formula in three dimensions become even more intractable.

Recall that the solution to Laplace's equation in one dimension is a straight line. Let us use Poisson's equation which is slightly more interesting:

$\nabla^2 u = \frac{\rho}{\epsilon_0}$

where $\rho(x)$ is the charge density (coulombs per cubic metre) and $\epsilon_0$ is the permittivity of free space (a constant). Both of these are given and therefore all we must do is find an approximation for $u(x)$.

We can rewrite this equation as $\frac{\partial^2}{\partial x^2}u(x) - \frac{\rho(x)}{\epsilon_0} = 0$

and thus, rather than solving directly for $u(x)$, we can define a function $V(x) = \frac{\partial^2}{\partial x^2}u(x) - \frac{\rho}{\epsilon_0}$ and thus our goal will be to find a function which satisfies $V(x) = 0$.

At this point, the reader is probably asking, "Yes, so what?". The point is that if $V(x) \equiv 0$ then it must also be true that

$\int_a^b V(x)\phi(x) dx = 0$

for any function $\phi(x)$ (we will call $\phi(x)$ a test function). What we will do is the following:

Partition the interval $[a, b]$ on which we wish to solve our boundary-value problem into $n$ sub-intervals such that $x_0 = a < x_1 < x_2 < \ldots < x_n = b$;
Find $n - 1$ test functions $\phi_1(x), \phi_2(x), \ldots, \phi_{n - 1}(x)$ such that we may use the formula $\int_a^b V(x)\phi_k(x) dx = 0$ to create a system of $n - 1$ equations and unknowns.

With a justification which will become apparent when we look into actually finding the linear systems, we will choose $\phi_k(x)$ to be that a continuous piecewise linear function, suggested by Galerkin, which is zero outside the interval $[x_{k - 1}, x_{k + 1}]$ and is 1 at $x_k$. Examples of such test functions are shown in Figure 3.

Figure 3. Various test functions $\phi_k(x)$.

Finding the equations of the sides of the triangle is straight-forward: the rising edge of the test function is $\frac{x - x_{k - 1}}{x_k - x_{k - 1}}$ while the falling edge is $\frac{x - x_{k + 1}}{x_k - x_{k + 1}}$. Therefore, our test functions are

$\phi_k(x) = \left \{ \begin{matrix} 0 & x \le x_{k - 1} \cr \frac{x - x_{k - 1}}{x_k - x_{k - 1}} & x_{k - 1} \le x \le x_k \cr \frac{x - x_{k + 1}}{x_k - x_{k + 1}} & x_k \le x \le x_{k + 1} \cr 0 & x \ge x_{k + 1} \end{matrix} \right .$.

(2)

(The operator ≤ is used everywhere to emphasize that the test function is continuous everywhere.)

Recall that our goal is to find points on the interior of the interval such that $u(x_k) \approx u_k$, as is shown in Figure 4.

Figure 4. The unknown interior points.

Note that we are not integrating $\phi_k(x)$ and $u(x)$, but rather we are integrating $\phi_k(x)$ and $\frac{d^2}{dx^2} u(x) - \frac{\rho(x)}{\epsilon_0}$ and we expect that integral to be zero.

The integral is

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x) - \frac{\rho(x)}{\epsilon_0} \right)\phi_k(x) dx = 0$.

Using the rules of integration, we have

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx - \int_{x_{k - 1}}^{x_{k + 1}} \frac{\rho(x)}{\epsilon_0} \phi_k(x) dx = 0$

where the second may be easily calculated because the density $\rho(x)$ is a given function. In the special case where $\rho(x)$ is constant (let $\frac{\rho}{\epsilon_0} = \rho_0$), the equation simplifies to

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx - \rho_0 \frac{x_{k + 1} - x_{k - 1}}{2} = 0$

(3)

and thus we can concentrate on the first integral without distraction. The biggest issue is approximating the second derivative! As discussed above, this requires reasonably well laid out points; especially in two and three dimensions. If we are to attempt to use tessellations as in Figure 2, this will be exceptionally difficult. Thus, we must attempt to get rid of the second derivative. We can do so by using integration by parts: $\int u'v = uv - \int uv'$ (found by integrating and manipulating the product rule $(uv)' = u'v + uv'$) only we will substitute $u$ with $u'$ and use the formulation $\int u''v = u'v - \int u'v'$; thereby getting rid of the second derivative.

The first integral may be rewritten using integration by parts as

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx = \left. \left(\frac{d}{dx} u(x)\right) \phi_k(x) \right |_{x = {x_k - 1}}^{x_{k + 1}} - \int_{x_{k - 1}}^{x_{k + 1}} \frac{d}{dx} u(x) \frac{d}{dx} \phi_k(x) dx $

(4)

You will immediately notice that we chose our test function $\phi_k(x)$ so that it is zero at both $x_{k - 1}$ and $x_{k + 1}$; therefore, the integral simplifies to

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx = -\int_{x_{k - 1}}^{x_{k + 1}} \frac{d}{dx} u(x) \frac{d}{dx} \phi_k(x) dx $.

The next simplification is that we can break the integral into two integrals:

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx = -\int_{x_{k - 1}}^{x_k} \frac{d}{dx} u(x) \frac{d}{dx} \phi_k(x) dx - \int_{x_k}^{x_{k + 1}} \frac{d}{dx} u(x) \frac{d}{dx} \phi_k(x) dx $

which therefore has the integral defined on each non-zero straight segment in Equation 2. At which point, we note that the derivatives of $\phi_k(x)$ in Equation 2 on each segment of integration are constant on each of the intervals. This can be seen in Figure 3 where the slopes of the two sides of $\phi_3(x)$ is easily calculated.

Figure 3. The slopes (red) of the two sides of the test function $\phi_3(x)$.

Substituting these two slopes into the integrands yields

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx = -\int_{x_{k - 1}}^{x_k} \frac{d}{dx} u(x) \frac{1}{x_k - x_{k - 1}} dx - \int_{x_k}^{x_{k + 1}} \frac{d}{dx} u(x) \frac{1}{x_k - x_{k + 1}} dx $.

The next step is to approximate the derivatives on each of the intervals. Once again, we can go to first-year calculus and we note that if we look at Figure 4, for example, the slope of the function $u(x)$ is well approximated on the interval $[x_2, x_3]$ by $\frac{u_3 - u_2}{x_3 - x_2}$, and the slope on the interval $[x_3, x_4]$ is well approximated by $\frac{u_4 - u_3}{x_4 - x_3}$.

Figure 4. Approximations (light blue) of the slope of $u(x)$ on the intervals $[x_2, x_3]$ and $[x_3, x_4]$.

Thus, we may substitute these two which we may do with the basic formula:

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx \approx -\int_{x_{k - 1}}^{x_k} \frac{u_k - u_{k - 1}}{x_k - x_{k - 1}} \frac{1}{x_k - x_{k - 1}} dx - \int_{x_k}^{x_{k + 1}} \frac{u_{k + 1} - u_k}{x_{k + 1} - x_k} \frac{1}{x_k - x_{k + 1}} dx $.

Now we note that both integrands are constants and we may use simple integration:

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx \approx -\frac{u_k - u_{k - 1}}{x_k - x_{k - 1}} \frac{1}{x_k - x_{k - 1}} (x_k - x_{k - 1}) - \frac{u_{k + 1} - u_k}{x_{k + 1} - x_k} \frac{1}{x_k - x_{k + 1}} (x_{k + 1} - x_k) $.

This how simplifies to

$\int_{x_{k - 1}}^{x_{k + 1}} \left( \frac{d^2}{dx^2} u(x)\right)\phi_k(x) dx \approx -\frac{u_k - u_{k - 1}}{x_k - x_{k - 1}} + \frac{u_{k + 1} - u_k}{x_{k + 1} - x_k} $.

Substituting this into Equation 3, we get

$-\frac{u_k - u_{k - 1}}{x_k - x_{k - 1}} + \frac{u_{k + 1} - u_k}{x_{k + 1} - x_k} - \rho_0 \frac{x_{k + 1} - x_{k - 1}}{2} = 0$.

We may rewrite this as

$\frac{1}{x_k - x_{k - 1}}u_{k - 1} - \left(\frac{1}{x_k - x_{k - 1}} + \frac{1}{x_{k + 1} - x_k}\right)u_k + \frac{1}{x_{k + 1} - x_k}u_{k + 1} = \rho_0 \frac{x_{k + 1} - x_{k - 1}}{2}$.

(5)

We now have a linear equation which we can now use.

Laboratory

First, we will look at the Matlab diff function and then we will continue with our function.

The `diff` Function

A function which you could easily use (and which would be very helpful) is the diff function. Given a vector v = [v(1), v(2), v(3), ..., v(n)], the command diff(v) returns an array of size $n - 1$ equal to [v(2) - v(1), v(3) - v(2), v(4) - v(3), ..., v(n) - v(n - 1)].

For example:

>> v = rand( 1, 5 )
v =
    0.1576    0.9706    0.9572    0.4854    0.8003
>> diff( v )
ans =
    0.8130   -0.0134   -0.4718    0.3149

If you call diff( v, 2 ), you get the result of diff( diff( v ) ), which is something which you will not need. However, you may need a vector of size $n - 2$ equal to [v(3) - v(1), v(4) - v(2), v(5) - v(3), ..., v(n) - v(n - 2)].

>> diff( v, 2 )
ans =
   -0.8264   -0.4584    0.7867
>> v(3:end) - v(1:(end - 2))
ans =
    0.7996   -0.4852   -0.1569

Program

Write the function

function u = elements1d( rho, xs, u_a, u_b )
    % Set up the matrix M
    % Set up the vector b

    % Update both the first and last entries of b
    b(1) =   b(1)   + ...;
    b(end) = b(end) + ...;
    
    u = [u_a; M \ b; u_b];
end

which implements our finite element method. In this laboratory, we will assume that the first argument rho is a constant. The argument xs is a sorted vector of values of the form [a, x(2), x(3), ..., x(n - 1), b] where $a < x_2 < x_3 < \cdots < x_{n - 1} < b$ and where $u(a) = u_a$ and where $u(b) = u_b$ are given. Thus, we are solving a system of length( xs ) - 2 equations in length( xs ) - 2 unknowns.

Questions

1. Apply your function to find the solution to

xs = linspace( 1, 6, 10 )';
u = elements1d( 0, xs, 3, 5 );
plot( xs, u );

Submit a plot of your solution.

2. Repeat Question 1 but use two random distributions of points.

xs1 = [1; (1 + 5*sort(rand( 8, 1 ))); 6];
xs2 = [1; (1 + 5*sort(rand( 8, 1 ))); 6];
u1 = elements1d( 0, xs1, 3, 5 );
u2 = elements1d( 0, xs2, 3, 5 );
plot( xs1, u1 );
hold on
plot( xs2, u2, 'r' );

Submit a plot of your solution.

3. Repeat Questions 1 and 2 but let the first argument, rho, be 1. Together with each of the two plots, include a plot of the actual solution in blue: $u(x) = 0.5x^2 - 3.1x + 5.6$. Use 100 points in your plot of the actual solution.

4. Repeat Questions 1 and 2 but let the first argument, rho, be -1. Together with each of the two plots, include a plot of the actual solution in blue: $u(x) = -0.5x^2 + 3.9x - 0.4$. Use 100 points in your plot of the actual solution.

5. Demonstrate that if we substitute Equation 1 into Poisson's equation, we get the same result we found in Equation 5. You may ask: "Why did we go through this whole effort to get something we could do half an hour ago?" The justification is that Equation 1 does not generalize to two and three dimensions. This technique of finite elements does as we will see in Laboratory 9.

NE 217 Advanced Calculus 2 for Nanotechnology Engineering