Numerical Analysis of Applied Partial Differential Equations

Bennett D. Ward

December 9, 2000

1  Introduction

Numerical analysis is the mathematical method that uses numerical approximations to obtain numerical answers to the problem. Numerical analysis also considers the accuracy of an approximation, and when the approximation is good enough. Numerical answers are useful because we use numbers to build our world, not with the exact analytical solution, such as e^p /Ö[27]

2  Numerical Approximations for Differential Calculus

Whereas the definition of the derivative is a difference quotient:

du dx = lim Dx ® 0  u(x+Dx) - u(x) Dx

(1)

Furthermore, there are approximations to definition of the derivative. Each has its own benefits and disadvantages.

du dx ® u(x+Dx) - u(x) Dx (2) du dx ® u(x) - u(x-Dx) Dx (3) du dx ® u(x+ 1 2 Dx) - u(x- 1 2 Dx) Dx (4)

where (2) is the forward difference equation, (3) is the backward difference equation, and (4) is the central difference. The forward difference is the difference equation that would be used on the left boundary, when the function proceeds to the right. This is most common when the function is set to exist in the first quadrant only. The backward difference would be used to evaluate functions at the opposite boundary, where the function would have no reference outside of this boundary. Therefore the difference may only be evaluated at the boundary and what we know about the function prior to that point.

The central difference is usable everywhere in between. The central difference requires information in front of, and behind the location being evaluated. Whereas the forward and backward differences are evaluated with the information on only one side of the location being evaluated.
For these reasons the central difference has a magnitude of error from the derivative of

1 24 (Dx)2 u(3)(x)

(5)

when Dx is small. The error of the forward difference is 1/2Dxu^¢¢(x). The magnitude of error for the backward difference ~ Dx The central difference is used as often as possible to obtain the most accurate results.

2.1  Higher Order Differential Calculus

d2u dx2 ® 1 Dx æ ç è du dx (x+Dx) - du dx (x) ö ÷ ø

(6)

when this is evaluated and expanded it becomes

d2u dx2 ® 1 Dx æ ç è u(x+2Dx)-u(x+Dx Dx - u(x+Dx)-u(x) Dx ö ÷ ø

(7)

This form is much more useful, and from it we can generalize the forward, central, and backward difference formulas for the second derivative. This also allows us the comfort of any combination of the difference formulas.

3  Double Ended Boundary Value Examples

d2u dx2 = f(x),   0 < x < 1 u(0) = 0,    u(1) = 1

Our central difference formula for the second derivative is

d2u dx2 ® u(x+Dx) - 2u(x) + u(x-Dx) (Dx)2

(8)

which we will use to solve our differential equation. If we specifically choose f(x) = -1, and say 4 partitions, Dx = 1/4 then we can build a system of equations which can easily be solved for the unknown values of the approximate functional value of u at the edges of the partitions. Our system of equations will then look like

d2u dx2 ® -1 = u(0) - 2u(1) + u(2) (1/4)2

(9)

d2u dx2 ® -1 = u(1) - 2u(2) + u(3) (1/4)2

(10)

d2u dx2 ® -1 = u(2) - 2u(3) + u(4) (1/4)2

(11)

replacing u(0) and u(1) with their known values, we can use any known method to solve for the three unknowns.

For the record they evaluate to u₁(x) = 11/32, u₂(x) = 5/8, and u₃(x) = 27/32.

If we change our differential equation to [(d²u)/(dx²)] -u = -2x and leave our other conditions the same, our system of equations is setup using the same method, where x is determined by Dx and the current interval being evaluated.

u0 - 2u1 + u2 (1/4)2 -u1 = -2(1*Dx) u1 - 2u2 + u3 (1/4)2 -u2 = -2(2*Dx) u2 - 2u3 + u4 (1/4)2 -u3 = -2(3*Dx)

then by replacing the known boundary conditions, (u₀ and u₄), we can solve by methods of linear equations.

If we had chosen to have more partitions, that is Dx is smallern ³ 4, then we would have one equation for each partition. Which would give a better representation of the actual graph. Each partition gives us the value of the function at the corresponding x location from information about the function on both sides of that function. Smaller Dx gives narrower partitions, which increases the accuracy of our approximation, but requires n-1 equations to evaluate, at best.
The definition of our difference formula determines the number of equations we must use, and how we must set it up. The more we know about our system, the less unknowns we have. In the two examples, we knew the functional value of u(x) at the endpoints, (u₀ and u₄), this allowed us to evaluate the remaining unknowns.

4  One sided BVP

ui+1-2ui+ui-1 (Dx)2 = f(xi) du dx = -1

for i = 1,2,3,...,n; where we have problems at

un+1-2un+un-1 (Dx)2

= f(xn)

(12) because we index past the right endpoint. But if we use the central difference, which we know, in place of the right endpoint, then we are covered and our approximate equation becomes

2un-1-2un -2Dx (Dx)2

= f(xn)

(13) With this, we can setup our system of equations for our number of unknowns and solve as before.

Differential equations of the form

d2u dx2 = f(x),    0 < x < 1 u(0)-gu¢(0) = c,     u(1) = 0

may be evaluated by the standard central difference formula that we have been using, Eq (8) to evaluate the homogeneous equation. For the non-homogeneous condition we will use (2), the forward difference formula, in place of u¢(0), our non-homogeneous condition is then

u0-g u1-u-1 2Dx = c

(14)

where u_-1 = 2Dx/g(c-u₀) which is then substituted into (14).

Not all BVP's have solutions, or the solutions do not fit near the boundary's. It is best to always check your answer with what you would expect. Sketching a graph from the numerical results is perhaps the simplest way to see how the numerical approximation fits with the anticipated solution, regardless of how the anticipated solution is derived.
For increased accuracy, a smaller Dx should be chosen. This will increase n, which increases the unknowns which ultimately results in better resolution of the graph plotted from the numerical results.

5  Numerical Methods for solving the Heat Equation

d2u dx2 = du dt

(15)

which is simply: second spatial derivative = first time derivative. Which will introduce two grid ``spacings''. The left hand side will remain with a Dx, but the right hand side will involve a Dt. We will see that this Dt has a maximum value for the solution to remain stable. The use of the forward difference in place of the time derivative is selected for two particular reasons. Time usually progresses forward, so we can approximate the next value from the data evaluated at the preceding grid location. This is done by solving for u(iDx+1,mDt+1) from the forward difference formula. If we replace (iDx+1,mDt+1) with (m+1), our heat equation in difference form is

Dt (Dx)2 [ui-1(m) -2ui(m)+ui+1(m)] = ui(m+1)-ui(m)

(16)

where i = 1,2,3,...,n-1 and m = 0,1,2,... if we solve this equation for (m+1) our equation then has the desired form of

ui(m+1) = Dt (Dx)2 ((ui-1(m)+ui+1(m)) )

(17)

which can then be setup in a system of equations to solve for the unknowns. These evaluations can then be used in the recursion to solve for the next set of unknowns. Whatever is known about a point on the grid may be placed into these equations. These solutions should always be checked graphically. Not all BVP's have solutions, and for a solution of this type of differential equation, the stability is dependent upon the grid size ratio ([(Dt)/((Dx)²)]).

If we choose our grid size ratio too big, then our solution will lose stability. By this, our values that are supposed to approximate our situation will fluctuate beyond ``tolerance'' and our approximation will be no good. On the other hand, if we arbitrarily choose a grid size too small, we may be overworking ourselves for no reason.

The hand waving method to choose the largest grid size ratio is as follows:

1. Write out the equations for each u_i(m) in terms of the previous time. such as we did in (17)
2. No coefficient is negative
3. The sum of the coefficients of u(m) is £ 1

Each differential equation, and BVP, will have its own maximum value for ([(Dt)/((Dx)²)]). However, this method should work for every problem. Numerical methods are needed most in situations with variable coefficients, inhomogeneities, or time dependent boundary conditions. These things must be considered when calculating the largest grid size ratio that will give a stable solution.

6  Use of Integrations

The limits of integration may be selected to fit the intervals of the desired difference equation. For example, to use a central difference the limits over each interval would be from the midpoint of the first interval, to the midpoint of the next interval. The intervals need not be of equal spacing.

This technique is also applicable when two materials are interfaced. Each material has its own solution to the temperature distribution. The distribution is continuous at the boundary, although the temperature may fall off very quickly.
If the temperature distributions follow the form [(¶)/(¶x)] ( k₁ [(¶u)/(¶x)]) = c₁[(¶u)/(¶t)] then the first partial may be eliminated through an integration wrt x. This leaves this problem in the form of k₁ [ [(u₂(m)-u₁(m))/(Dx₁)] - [(u₁(m)-u₀(m))/(Dx₁)] ] = c₁Dx₁[(u₁(m+1)-u₁(m))/(Dt)] where the subscript refers to the material.

This simplifies the form so that replacement equations may be used with less confusion. If the solution could be solved by a double integration then the difference equations could be avoided all together. However, the fact that we are required to use these difference methods suggests that an analytic solution will not be found by sequential integrations. Vice versa, the absence of analytical solutions is perhaps what suggested the use of the difference approximations.

From this point, the difference approximations may be used to supply the replacement equations, which in turn are evaluated to reveal the unknowns. We are already familiar with those methods.

7  Using what I've learned

d2u dx2 = du dt ,     u(0,t) = u(1,t) = t , u(x,0) = 0)

Using the methods illustrated in (16) and (17), our replacement equation satisfies the second condition for finding the grid size ratio. Having this condition satisfied, I will solve the equation for ([(Dt)/((Dx)²)]). This will ensure that the third condition is satisfied. I must note that I have considered the steps of my technique to find the grid size ratio as conditions. These steps must be conditions until I learn of another method to solve for this value.

assuming that we have pre-selected our Dx (from Dx = 1/n) then our Dt will evaluate to the largest stable time step. This ratio is calculated from an intermediate form of (16). u_i(m+1) = ru_i-1(m) +(1-2r)u_i(m) +ru_i+1(m) where r has been substituted for ([(Dt)/((Dx)²)]). Evaluating this intermediate form at n = 4 for r gives r = 1/2. This primarily happens due to all the coefficients =1 when represented in the form of (17). Our Dt then may be evaluated to 1/32

Our replacement equations will be

u1(m+1) = r(u0(m)+u2(m) ) u2(m+1) = r(u1(m)+u2(m) ) u3(m+1) = r(u2(m)+u4(m) ) (18)

This will represent the distribution over the interval, in the x direction, for a given t_m which is represented by mDt. Therefore we will increment m = 1,2,... because u_i(0) = 0.

At m ¹ 0, u₀(m) = u₄(m) = t This considers the interval endpoints of x at any m ¹ 0, where t = mDt. From this we may setup a spreadsheet, or linear algebra method to solve for everything in between.

m	0	1	2	3	4
0	0	0	0	0	0
1	0.03125	0	0	0	0.03125
2	0.0625	0.015625	0	0.015625	0.0625
3	0.09375	0.03125	0.015625	0.03125	0.09375
4	0.125	0.0546875	0.03125	0.0546875	0.125
5	0.15625	0.078125	0.0546875	0.078125	0.15625

Back to Mathlinks

Back to the Front