Approximate and Numerical Methods

From Scholarpedia

Curator: Dr. Andrei D. Polyanin, Institute for Problems in Mechanics, Moscow, Russia
Curator: Dr. William E. Schiesser, Lehigh University, USA
Curator: Dr. Alexei I. Zhurov, Cardiff University, UK, and Institute for Problems in Mechanics, Moscow, Russia.

First-Order Partial Differential Equations

Second-Order Partial Differential Equations

Higher-Order Partial Differential Equations

Approximate and Numerical Methods

The preceding discussion pertains to the exact or analytical solution of PDEs. For example, in the case of a heat equation or a wave equation, an exact solution would be a function which, when substituted into the respective equation would satisfy it identically along with all of the associated initial and boundary conditions.

Although analytical solutions are exact, they also may not be available, simply because we do not know how to derive such solutions. This could be because the PDE system has too many PDEs, or they are too complicated, e.g., nonlinear, or both, to be amenable to analytical solution. In this case, we may have to resort to an approximate solution. That is, we seek an analytical or numerical approximation to the exact solution.

Perturbation methods are an important subset of approximate analytical methods. They may be applied if the problem involves small (or large) parameters, which are used for constructing solutions in the form of asymptotic expansions. For books on perturbation methods, see Google Book Search. These and other methods for PDEs are also outlined in Zwillinger (1997).

Unlike exact and approximate analytical methods, methods to compute numerical PDE solutions are in principle not limited by the number or complexity of the PDEs. This generality combined with the availability of high performance computers makes the calculation of numerical solutions feasible for a broad spectrum of PDEs (such as the Navier–Stokes equations) that are beyond analysis by analytical methods. The development and implementation (as computer codes) of numerical methods or algorithms for PDE systems is a very active area of research. Here we indicate in the external links just two readily available links to Scholarpedia.

We now consider the numerical solution of a parabolic PDE, a hyperbolic PDE, and an elliptic PDE.

Parabolic PDE

Analytical solutions to a parabolic PDE (heat equation) are given here. But we will proceed with a numerical solution and use one of these analytical solutions to evaluate the numerical solution.

We can consider the numerical solution to the heat equation

(1)

as a two-step process:

1. Numerical approximation of the derivative . At this point, we will have a semi-discretization of Eq. (1).
2. Numerical approximation of the derivative . At this point, we will have a full discretization of Eq. (1).

In order to implement these two steps, we require a grid in and a grid in . For the grid in , we denote a position along the grid with the index . Then we can consider the Taylor series expansion of the numerical solution at grid point

(2)

and

(3)

If we consider a uniform grid (a grid with uniform spacing , addition of Eqs. (2) and (3) gives (note the cancellation of the first and third derivative terms since )

(4)

where denotes a term proportional to or of order ; this term can be considered a truncation error resulting from truncating the Taylor series of Eqs. (2) and (3) beyond the term. Then Eq. (4) gives for the second derivative

(5)

Equation (5) is a second order (because of the principal error or truncation error ) finite difference approximation of .

If Eq. (5) is substituted in Eq. (1) (to replace the derivative ), a system of ODEs results

(6)

(we have added a multiplying constant to the right-hand side of Eq. (1), generally termed a thermal diffusivity if in Eq. (1) is temperature and a mass diffusivity if is concentration; has the MKS units m²/s as expected from a consideration of Eq. (1) with in metres and in seconds).

Note that the independent variable does not appear explicitly in Eqs. (6) and that the only independent variable is (so that they are ODEs). is the number of points in the grid ( is termed a boundary value variable since the terminal grid points at and typically refer to the boundaries of a physical system). Thus Eq. (1) is partly discretized (in ) and therefore Eqs. (6) are referred to as a semi-discretization.

To compute a solution to Eq. (1), we could apply an established initial-value integrator in . This is the essence of the method of lines (MOL). Alternatively, we could now discretize Eqs. (6). For example, if we apply Eq. (2) on a grid in with an index

(7)

If the grid in has a uniform spacing and if truncation after the first derivative term is applied,

(8)

Equation (8), the classical Euler's method, can be used to step along the solution of Eq. (1) from point to (at a grid point in ).

Application of Eq. (8) to Eq. (6) gives the fully discretized approximation of Eq. (1)

(9)

In Eq. (7) we do not specify the total number of grid points in (as we did with the grid in ); is an initial value variable since it is typically time, and is defined over the seminfinite interval .

Note that Eq. (9) explicitly gives the solution at the advanced point in (at ) and therefore it is an explicit finite difference approximation of Eq. (1). We can now consider using Eq. (9) to step forward from an initial condition (IC) required by Eq. (1). Here we take as the initial condition

(10)

where are constants to be specified. The finite difference form of Eq. (10) is

(11)

(note that at ).

We must also specify two boundary conditions (BCs) for Eq. (1) (since it is second order in ). We will use the Dirichlet BC at

(12)

for which the finite difference form is (note that at )

(13)

We use the Neumann BC at

(14)

for which the finite difference form is (note that at )

(15)

where is a fictitious value that is outside the interval ; it can be used to eliminate in Eq. (9) for .

Equations (9), (11), (13) and (15) constitute the full system of equations for the calculation of the numerical solution to Eq. (1). Note that we have replaced the original PDE, Eq. (1), with a set of approximating algebraic equations (Eqs. (9), (11), (13) and (15)) which can easily be programmed for a computer.

Also, an analytical solution to Eq. (1) (see particular solutions to the heat equation) can be used to evaluate the numerical solution

(16)

Equation (16) can be stated in the alternative form

(17)

which corresponds to a traveling wave solution since and appear in the combination .

A short MATLAB program is listed in Appendix 1 based on Eqs. (9), (11), (13) and (15). Representative output from this program that compares the numerical solution from Eqs. (9), (11), (13) and (15) with the analytical solution, Eq. (17), indicates that the two solutions are in agreement to five figures, as reflected in Table 1.

Table 1: Comparison of the numerical and analytical solutions at produced by the program in Appendix 1; w, numerical solution from Eqs. (9), (11), (13) and (15); wa, analytical solution of Eq. (17)

t = 0.05   w(x=xl/2,t) = 1.6376   wa(x=xl/2,t) = 1.6376

t = 0.10   w(x=xl/2,t) = 1.6703   wa(x=xl/2,t) = 1.6703

t = 0.15   w(x=xl/2,t) = 1.7047   wa(x=xl/2,t) = 1.7047

t = 0.20   w(x=xl/2,t) = 1.7408   wa(x=xl/2,t) = 1.7408

t = 0.25   w(x=xl/2,t) = 1.7788   wa(x=xl/2,t) = 1.7788

t = 0.30   w(x=xl/2,t) = 1.8187   wa(x=xl/2,t) = 1.8187

t = 0.35   w(x=xl/2,t) = 1.8607   wa(x=xl/2,t) = 1.8607

t = 0.40   w(x=xl/2,t) = 1.9048   wa(x=xl/2,t) = 1.9048

t = 0.45   w(x=xl/2,t) = 1.9512   wa(x=xl/2,t) = 1.9512

t = 0.50   w(x=xl/2,t) = 2.0000   wa(x=xl/2,t) = 2.0000

The parameters that produced the numerical output in Table 1 are listed in Table 2.

Table 2: Numerical values of parameters used in the program of Appendix 1 that produced the output of Table 1

Parameter	Description	Value
	diffusivity in Eq. (1)	1
, ,	constants in Eqs. (10)-(17)	1, 1, 1
,	number of grid points and length in	21, 1
,	grid spacing in , final value of	0.001, 0.5

Additional parameters follow from the values in Table 2. Thus, the grid spacing in is . The number of steps in taken along the solution is .

Some of the parameters, particularly and , were determined by trial and error to achieve a numerical solution of acceptable accuracy (e.g., five significant figures). We can note two additional points about these values:

• Acceptable values of and could be determined by observing the errors in the numerical solution through comparison with the exact solution (Eq. (17)) as illustrated in Table 1. However, in most PDE applications, an analytical solution is not available for assessing the accuracy of the numerical solution, and in fact, the motivation for using a numerical method is generally to produce a solution when an analytical solution is not available. In this case (no analytical solution), a useful procedure for estimating the numerical accuracy is to compute solutions for two different values of and compare the numerical values. If the two solutions do not agree to an acceptable level, a third solution is computed with a still larger (smaller grid spacing in ) and again the solutions are compared. Eventually, if spatial convergence is achieved, successive solutions will agree to the required accuracy. In this case, the accuracy of the numerical solution is inferred, but the exact error is not computed (or even known) when an analytical solution is not available. The same reasoning can be applied for temporal convergence with respect to , i.e., is reduced until the successive solutions agree to a specified level.

• The preceding discussion was directed to achieving acceptable accuracy. Additionally, the values of and were selected to achieve a stable solution. The criterion for stability in the case of Eq. (9) (for parabolic PDE (1)) is . For the solution of Table 1, . If the critical value of the dimensionless parameter had been exceeded, the numerical solution would have become unstable (as manifest in numerical values of ever increasing magnitude). This conclusion can easily be confirmed by running the program of Appendix 1 for a choice of and for which . The stability constraint is a distinctive feature of the explicit finite difference approximation of Eq. (9). Thus, as is reduced ( increased) to achieve better accuracy in the numerical solution, must also be reduced to maintain stability ().

Finally, some descrpitive comments about the details of the program in Appendix 1 are given immediately after the program listing.

Hyperbolic PDE

We now consider a numerical solution to the one-dimensional hyperbolic wave equation

(18)

We again have an analytical solution to evaluate the numerical solution. First, we include a velocity, , in the equation:

(19)

Note that has the MKS units of m/s as expected and as inferred from Eq. (19) (note the units of the derivatives in and ).

Since Eq. (19) is second order in and , it requires two ICs and two BCs. We will take these as:

(20)

(21)

(22)

(23)

Equation (20) indicates the IC is a Gaussian pulse with the positive constant to be specified. Equation (21) indicates starts with zero "velocity". Equations (22)–(23) indicate that the solution does not depart from the initial value of zero specified by IC (20)–(21). In other words, is chosen large enough that the IC is effectively zero at and remains at this value for subsequent .

An important difference between the parabolic problem of Eq. (1) and the hyperbolic problem of Eq. (19) is that the former is first order in while the latter is second order in . Therefore, in order to develop a numerical method for Eq. (18) (or Eq. (19)), we need an algorithm that can accommodate second derivatives in . While such algorithms do exist, they generally are not required. Rather, we can express a PDE second order in as two PDEs first order in . For example, Eq. (19) can be written as

(24)

(25)

Equations (24)–(25) are first order in , and therefore an integration algorithm for first order equations, such as the Euler method of Eq. (8), can be used to move and forward in . Thus, the fully discretized form of Eqs. (24)–(25) can be written as

(26)

(27)

ICs (20)–(21) become (with corresponding to )

(28)

(29)

For BCs (22)–(23), the infinite interval must be replaced by a finite one (since computers can accommodate only finite numbers) where is selected so that it is effectively infinite; that is, the solution does not depart from IC (20) at for . The value of and the corresponding number of grid points in , , are specified subsequently.

The finite difference approximations of BCs (22)–(23) are

(30)

(31)

Equations (26), (27), (28), (29), (30), (31) constitute the complete finite difference approximation of Eqs. (24), (25), (20), (21), (22), (23). A small MATLAB program for this approximation is given in Appendix 2, including some descriptive comments immediately after the program.

The general analytical solution to Eq. (19) can be written as (see also here)

(32)

Let us take in the form of the Gaussian pulse of Eq. (20), i.e.,

(33)

Figure 1: Comparison of the numerical and analytical solutions for produced by the program in Appendix 2; w, numerical solution from Eqs. (26), (27), (28), (29), (30), (31); wa, analytical solution of Eq. (33)

The plotted output from the program is given in Figure 1 and includes both the numerical solution of Eqs. (26), (27), (28), (29), (30), (31) and the analytical solution of Eq. (33).

We can note the following points about Figure 1:

• The initial Gaussian pulse at (centered at with unit maxiumum value) splits into two pulses traveling left and right with velocity and maximum value of 0.5 according to Eq. (33).

• The pulses traveling left are centered at corresponding to since . The pulses traveling right are centered at corresponding to . These properties are characteristic of the traveling wave functions of Eq. (33) with arguments and , respectively. In fact, the use of the word characteristic is particularly appropriate since the relations and , where is a constant, are termed the characteristics of Eq. (19). Note that at the points along the solution given by these characteristics, the solution has a constant value. For example, for the peak values in Figure 1, and the solution is constant at the peak value 0.5 for .

• The solution remains at zero for so that the interval is equivalent to the infinite interval .

The parameters that produced the numerical output in Figure 1 are listed in Table 3.

Table 3: Numerical values of parameters used in the program of Appendix 2 that produced the output of Figure 1

Parameter	Description	Value
	velocity in Eq. (19)	1
	constant in Eqs. (20), (28), (33)	0.05
,	number of grid points and half length in	201, 50
,	grid spacing in , final value of	0.0025, 30

Additional parameters follow from the values in Table 3. Thus, the grid spacing in is . The number of steps in taken along the solution is , which is large, but was selected to achieve good accuracy in the numerical solution.

To explore the accuracy of the numerical solution, we can consider the peak values of in Figure 1. This is a stringent test of the numerical solution since the curvature of the solution is greatest at these peaks. The results are summarized in Table 4.

Table 4: Numerical values of at the peak values displayed in Figure 1

		Peak values
0	0	1
10	-10, 10	0.5006, 0.5006
20	-20, 20	0.5011, 0.5011
30	-30, 30	0.5015, 0.5015

Of course, the peak analytical values given by Eq. (33) are , . Table 4 indicates these peak values were attained within 0.5015 even for the largest value of so that errors did not accumulate excessively as the solution progressed through . These errors could be reduced by increasing the number of grid points in above . These errors could also presumably be reduced by using a more accurate (higher order) finite difference equation than the second order approximation of Eq. (5); to this end, MATLAB routines for second, fourth, sixth, eighth and tenth finite difference approximations are given in Hamdi, et al.

This explicit finite difference numerical solution also has a stability limit (like the preceding parabolic problem). In this case

(34)

Stability constraint (34) is the Courant–Friedrichs–Lewy (or CFL) condition. For the present numerical solution, so the CFL condition is easily satisfied. In other words, the parameters of Table 3 were chosen primarily for accuracy and not stability.

Next, we consider an elliptic problem.

Elliptic PDE

The elliptic PDE (Laplace's equation)

(35)

has two boundary value independent variables, and , and no initial value variable. Thus, since the preceding numerical methods required an integration with respect to an initial value variable (), and was accomplished by Euler's method, Eq. (8), we cannot develop a numerical solution for Eq. (35) directly using these methods. Rather, we will convert Eq. (35) from an elliptic problem to an associated parabolic problem by adding a derivative in an initial value variable. Thus, the PDE we will now consider is

(36)

The idea then is to integrate Eq. (36) forward in until the numerical solution approaches the condition so that Eq. (36) then reverts back to Eq. (35), i.e., the solution under this condition is for Eq. (35) as required.

Equation (36) is first order in , and second order in and . It therefore requires one IC and two BCs (for and ). For the IC, since has been added to the problem only to provide a solution to Eq. (35) through Eq. (36), the choice of an IC is completely arbitrary (it is not part of the original problem). For the present analysis, we will use

(37)

where is a constant to be selected (logically, it should be in the neighborhood of the expected solution to Eq. (35), but it's precise value is not critical to the success of the numerical method).

For the two BCs in , we will use homogeneous (zero) Dirichlet BCs

(38)

(39)

where is the upper boundary value of .

For the first BC in , we will use a homogeneous Neumann BC

(40)

For the second BC in , we will use a nonhomogeneous Neumann BC

(41)

where is the upper boundary value of . Note that in BC (41) is a function of .

The analytical solution to Eqs. (35) (note, not Eq. (36)), (38)–(41) is a special case of one of the solutions stated previously

(42)

The analytical solution of Eq. (42) will be used to evaluate the numerical solution.

To develop a numerical solution to Eq. (36), we first replace all of the derivatives with finite differences in analogy with the preceding numerical solutions. The positions in , and will be denoted with indices , and , respecively. The corresponding increments are , and . Application of these ideas to Eq. (36) gives the finite difference approximation

(43)

BCs (38)–(39) in the difference notation are

(44)

(45)

BC (40) in difference notation is

(46)

where is a fictitious value that can be used in Eq. (43) for .

BC (41) in difference notation is

(47)

where is a fictitious value that can be used in Eq. (43) for .

A short MATLAB program for Eqs. (43), (44), (45), (46) and (47) is listed in Appendix 3, followed by some explanatory comments. The output from this program is listed in Table 5.

Table 5: Numerical values of from the program in Appendix 3

0.10	1.6254
0.20	2.1785
0.30	2.3883
0.40	2.4664
0.50	2.4956
0.60	2.5064
0.70	2.5105
0.80	2.5120
0.90	2.5125
1.00	2.5127

The midpoint values are listed in Table 5. The convergence of the solution of Eq. (36) to that of Eq. (35) is apparent, e.g., at , while the value from the analytical solution, Eq. (42), is ; the difference is . Of course, the agreement will not be perfect because of the approximations used in Eqs. (43), (44), (45), (46) and (47). The parameters that produced the numerical output in Table 5 are listed in Table 6.

Table 6: Numerical values of parameters used in the program of Appendix 3 that produced the output of Table 5

Parameter	Description	Value
,	number of grid points and half length in	21, 1
,	number of grid points and half length in	21, 1
,	grid spacing in , final value of	0.0005, 1

Additional parameters follow from the values in Table 6. Thus, the grid spacing in is . Similarly the spacing in is . The spacing in and could be different if the variation of the solution in one direction is substantially greater than in the other; in other words, some tuning of the 2D grid in and could be required. The number of steps in taken along the solution is , which is large, but was selected to achieve good accuracy in the numerical solution.

The stability constraint for the 2D problem of Eq. (36) (with a constant multiplying the derivatives in and ) is

For the preceding solution, which meets the stability constraint.

The procedure of adding a derivative in to elliptic Eq. (35), then solving the resulting parabolic problem, Eq. (36), is termed the method of false transients or the method of pseudo transients to suggest that the parabolic problem appears to have a transient phase that is not part of the original elliptic problem. The actual path that the parabolic problem takes to the solution of the elliptic problem is not relevant so long as the parabolic solution converges to the elliptic solution. In the present case, in IC (37) was taken as 1. Some other trial values indicated that this initial value is not critical (but it should be as close to the final value, for example 2.5092, as knowledge about the final value can provide).

Also, can be considered a parameter in the solution of Eq. (36). This parameterization is an example of continuation in which the solution is continued from the given (assumed) starting value of Eq. (37) to the final value of 2.5127. The concept of continuation has been applied in many forms (and not just through the addition of a derivative as in Eq. (36)).

In general, the errors in the numerical solution of PDEs can result from the limited accuracy of all of the approximations used in the calculation. For example, the 0.14% error in the preceding solution can result from the discretzation errors in , and in Eqs. (43), (44), (45), (46) and (47). In formulating a numerical method or algorithm for the solution of a PDE problem, it is necessary to balance the discretization errors so that one source of error does not dominate, and generally degrade, the numerical solution. This might, for example, require the choice of balanced values for , and . The effect of the numbers of discretization points generally can be inferred by varying these numbers and observing the effect on the numerical solution.

Remarks

Thus, control of approximation errors is central to the calculation of a numerical solution of acceptable accuracy. In the preceding examples, this control of errors can be accomplished in three ways:

1. The number of grid points can be varied. Since in the numerical analysis literature, the grid spacing is often given the symbol (as in Eq. (8)), systematic variation of the grid spacing to investigate accuracy is usually termed h refinement.
2. The order of the approximation can be varied. For example, Euler's method of Eq. (8) is (second order correct) for one step or (first order correct) for a series of steps over a complete solution. In general, if the approximation is of order , it is termed th order. The symbol is frequently used in the numerical analysis literature, and varying the order to investigate accuracy is referred to as p refinement.
3. The discretization errors can possibly be estimated and where the error is considered too large, the numerical algorithm can automatically insert grid points. Beyond that, the grid spacing does not have to be uniform, e.g., in Eq. (8) does not have to remain constant throughout the numerical solution, and the numerical algorithm can automatically vary the spacing to concentrate the grid points where they are needed to achieve acceptable accuracy. This form of refinement is termed usually r refinement; we mention r refinement only to initiate possible interest in more advanced numerical methods.

The three preceding numerical solutions were developed using basic finite differences such as in Eqs. (5) and (8). However, many approaches to approximating derivatives in PDEs have been developed and used. Among these are finite elements, finite volumes, weighted residuals, e.g., collocation, Galerkin and spectral methods. Each of these methods has advantages and disadvantages, often according to the characteristics of the problem of interest (starting with the parabolic, hyperbolic and elliptic geometric classifications). Thus, an extensive literature for the numerical solution of PDEs is available, and we have only presented here a few basic concepts and examples.

The principal advantage of numerical methods applied to PDEs is that, in principle, PDEs of any number and complexity can be solved which is particularly useful when analytical solutions are not available. To illustrate how this might be done, the heat conduction equation with a nonlinear source term can easily be programmed by extending Eq. (9) to

(48)

where could be a MATLAB function to compute the nonlinearity in . Alternatively, the nonlinearity could be programmed directly in Eq. (48). For example, for a second order reaction,

(49)

Note in particular how easily the nonlinearity is programmed.

As another example, a solution to the Burgers equation could be computed by extending Eq. (9) (with ) to

(50)

where we have used the second order central finite difference approximation

Again, the nonlinear term in the Burgers equation is easily programmed as

Finally, to conclude this introductory discussion of the numerical solution of PDEs, we used in the various examples integration with respect to an initial value variable, , by using the Euler method of Eq. (8). While Euler's method is general with respect to the form of the initial value integration, it does have two important limitations:

• The accuracy is first order () and is therefore exact only for functions that vary linearly in ; for higher order variations in , it is only approximate and generally requires a small to achieve acceptable accuracy. One way around this accuracy limitation is to use a higher order integration algorithm in such as the Runge-Kutta method (e.g., see the article by Shampine et al).

• Stability limits the step , such as in the stability constraint for Eq. (9), . Such stability constraints for explicit integration methods such as Eq. (8) can be circumvented by using an implicit integration method, often termed a stiff method. We will not pursue various types of initial value integration methods here, but additional information is available (Shampine et al).

Thus, the Euler method is limited by both accuracy and stability. To circumvent these constraints, initial value integrators are used which are higher order (for improved accuracy) and stable (to provide a larger stable step that may actually be unlimited, depending on the integrator and the application). As might be expected, such higher order methods are more complicated than the Euler method, but fortunately, they have been programmed in library routines that can easily be called and used. As an example, the MATLAB integrator ode15s has the following features:

• Good stability (it is a stiff integrator). The enhanced stability is achieved through implicit integration that requires the solution of algebraic equations at each step along the solution (each step of length ). An option for ode15s in the solution of the algebraic equations is the use of sparse matrix methods that are very efficient for large algebraic systems.

• Variable step and order, so that it automatically performs and refinement as the solution evolves.

• Library routines such as ode15s are written by experts who have included features that make them robust and reliable. Further, these routines have been used extensively and are therefore thoroughly tested.

The use of library routines for initial value integration is the basis for much of the work in the numerical method of lines solution of PDEs. This approach has been applied to a broad spectrum of PDE problems in 1D, 2D and 3D, including all of the major classes of PDEs (e.g., parabolic, hyperbolic and elliptic) in various orthogonal coordinate systems (e.g., Cartesian, cylindrical and spherical). The use of quality library routines provides an important step in the timely development of computer codes for new applications of PDEs whereby the analyst can take advantage of the work of experts, which is generally much more efficient and reliable than developing codes starting with just a general programming language.

Appendix 1. MATLAB program for a parabolic PDE

The MATLAB program that produced the results described in the Parabolic PDE for equation (1) follows.

 clear all; clc
 D=1; mu=1; A=1; B=1;
 nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11;
 nt=50; nout=10; h=0.001; t=0;
 for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
 for i1=1:nout
   for i2=1:nt
     for i=1:nx
       if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;
       elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
                    wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;
       else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
     end
     w=w+wt*h; t=t+h;
   end
   fprintf('t = %4.2f   w(x=xl/2,t) = %6.4f   wa(x=xl/2,t) = %6.4f\n', ...
           t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
 end

Listing 1: MATLAB program for the numerical solution of Eqs. (1), (10), (12), and (14)

We can note the following points about this program:

• Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (1) are defined numerically as in Table 2. Note that a grid in is set up as .

 clear all; clc
 D=1; mu=1; A=1; B=1;
 nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11;
 nt=50; nout=10; h=0.001; t=0;

• IC (11) is then implemented in a for loop (note that t=0).

 for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end

• Three nested for loops step Eq. (9) through and

 for i1=1:nout
   for i2=1:nt
     for i=1:nx

Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 1. Within each output interval (0.05), 50 steps in are taken by the loop in i2 (in this way, the 500 steps in produce outputs at only 10 points as displayed in Table 1). Within the loop in i, spans the interval .

• BC (13) is programmed (at i=1) as

       if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;

Also, since is specified by BC (12), the derivative of this boundary value is set to zero so that the integration in at does not move from its prescribed value, i.e., wt(1)=0;).

• BC (15) is used to calculate the fictitious point wf which is then used in the finite difference calculation of Eq. (6).

       elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
                    wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;

• For the interior points in (), Eq. (6) is programmed.

       else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
     end

The second end terminates the inner loop in i, so that all of the values of are handled.

• Euler's method, Eq. (8), is then used to advance the solution through nt=50 steps of length through the intermediate for loop in i2.

     w=w+wt*h; t=t+h;
   end

Note the use of the matrix facility of MATLAB (w and wt are vectors of length nx = 21). The end completes the intermediate loop in i2.

• The numerical solution is displayed to produce the output in Table 1, including the analytical solution of Eq. (16)

   fprintf('t = %4.2f   w(x=xl/2,t) = %6.4f   wa(x=xl/2,t) = %6.4f\n', ...
           t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
 end

The end terminates the outer loop in i1.

Appendix 2. MATLAB program for a hyperbolic PDE

The MATLAB program that produced the results described in the Hyperbolic PDE section for equation (18) (and Eq. (19)) follows.

 clear all; clc
 c=1;
 nx=201; xl=50; dx=0.5; x=[-xl:dx:xl];
 nt=4000; nout=3; h=0.0025; t=0;
 for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
             wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end
 for i1=1:nout
   for i2=1:nt
     for i=1:nx
       if(i==1)     w( 1)=0; wt( 1)=0; wtt( 1)=0;
       elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;
       else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
     end
     w=w+wt*h; wt=wt+wtt*h; t=t+h;
   end
 for i=1:nx
   wplot (i1+1,i)=w(i);
   waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
   fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
 end
 end
 plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
 title('w(x,t), t=0,10,20,30, solid - num, o - anal')

Listing 2: MATLAB program for the numerical solution of Eqs. (19), (20)–(23)

We can note the following points about this program:

• Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (18) (or Eq. (19)) are defined numerically as in Table 3. Note that a grid in is set up as , a total of 201 points. The choice of the upper and lower limits for , was made so that domain in is essentially infinite (it accommodates BCs (22)–(23), (30)–(31) as discussed subsequently).

 clear all; clc
 c=1;
 nx=201; xl=50; dx=0.5; x=[-xl:dx:xl];
 nt=4000; nout=3; h=0.0025; t=0;

• ICs (20)–(21), or (28)–(29), are then implemented in a for loop (note that t=0).

 for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
             wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end

Both and are programmed (as w and wt) which is required since Eq. (19) is second order in (as explained subsequently). Also, the initial value of , both numerical and analytical, is stored for subsequent plotting (so that the plot will include IC (20)).

• Three nested for loops step Eqs. (26) and (27) through and

 for i1=1:nout
   for i2=1:nt
     for i=1:nx

Specifically, the outer loop (in it) gives the nout = 3 outputs in Table 3. Within each output interval (1), 4000 steps in are taken by the loop in i2 (in this way, the 4000 steps in produce outputs at only 3 points as displayed in Figure 1). Within the loop in i, spans the interval .

• BCs (22)–(23), or (30)–(31), are programmed (at i=1 and i=nx corresponding to ) as

       if(i==1)     w( 1)=0; wt( 1)=0; wtt( 1)=0;
       elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;

Also, since is specified by BCs (22)–(23), the first and second derivatives in of these boundary values are set to zero so that the integration in does not move from its zero values.

• For the interior points in (), Eq. (25) (discretized in ) is programmed.

       else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
     end

The second end terminates the inner loop in i, so that all of the values of are handled.

• Euler's method (Eq. (8)) is then used to advance the solution through nt=4000 steps of length through the intermediate for loop in i2 according to Eqs. (26) and (27).

    w=w+wt*h; wt=wt+wtt*h; t=t+h;
   end

Note the use of the matrix facility of MATLAB (w, wt and wtt are vectors of length nx = 201). The end completes the intermediate loop in i2.

• The numerical and analytical solutions are then stored at (three passes through the outer loop in i1) followed by plotting (of Figure 1).

 for i=1:nx
   wplot (i1+1,i)=w(i);
   waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
   fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
 end
 end
 plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
 title('w(x,t), t=0,10,20,30, solid - num, o - anal')

The final end terminates the outer loop in i1. Also, as a word of caution, the fprintf statement displays the complete numerical solution at so that the numerical solution can be examined in detail (the peak values reported in Table 4 were taken from this output); if the program in Listing 2 is executed, the fprintf statement could be converted to a comment to reduce the numerical output.

Appendix 3. MATLAB program for an elliptic PDE

The MATLAB program that produced the results described in the Elliptic PDE section for equation (35) follows.

 clear all; clc
 nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl;
 ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11;
 nt=200; nout=10; h=0.0005; t=0;
 w=ones(nx,ny);
 for i1=1:nout
   for i2=1:nt
     for i=1:nx
       for j=1:ny
         if(i== 1)      w( 1,j)=0; wt( 1,j)=0;
           elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;
           elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                                +2*(w(i,2)-w(i,1))/dy^2;
           elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
                        wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
                                +(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
           else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                       +(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
         end
       end
     end
     w=w+wt*h;
     t=t+h;
   end
   fprintf('t = %5.2f   w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
 end
 fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));

Listing 3: MATLAB program for the numerical solution of Eqs. (36)–(41)

We can note the following points about this program:

• Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (36) are defined numerically as in Table 6. Note grids in and are set up as , .

 clear all; clc
 nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl;
 ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11;
 nt=200; nout=10; h=0.0005; t=0;

• IC (37) is then implemented as

 w=ones(nx,ny);

in Eq. (37) is programmed with the MATLAB utility ones over the entire grid in and , a total of nxny points.

• Four nested for loops step Eq. (43) through , and

 for i1=1:nout
   for i2=1:nt
     for i=1:nx
       for j=1:ny

Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 5. Within each output interval (0.0005), 200 steps in are taken by the loop in i2 (in this way, the 2000 steps in produce outputs at only 10 points as displayed in Table 5). Within the loop in i, spans the interval and within the loop in j, spans the interval .

• BC (38) is programmed (at i=1 for j=1 to j=ny) as

         if(i== 1)      w( 1,j)=0; wt( 1,j)=0;

Also, since is specified by BC (38), the derivative of this boundary value is set to zero so that the integration in at does not move from its prescribed value, i.e., wt(1,j)=0;).

• BC (39) is programmed (at i=nx for j=1 to j=ny) as

           elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;

• Equation (43) (with subsequent application of Euler's method) is programmed using BC (40) (at i=1 to i=nx for j=1) as

          elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                               +2*(w(i,2)-w(i,1))/dy^2;

Note that we have used a finite difference approximation of Eq. (40), w(i,0) = w(i,2), to eliminate the fictitious value w(i,0) in Eq. (43).

• Equation (43) (with subsequent application of Euler's method) is programmed using BC (41) (at i=1 to i=nx for j=ny) as

          elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
                        wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
                                +(wf-2*w(i,ny)+w(i,ny-1))/dy^2;

Note that we have used a finite difference approximation of Eq. (41) to eliminate the fictitious value w(i,ny+1)=wf in Eq. (43).

• For the interior points in () and (), Eq. (43) (with subsequent application of Euler's method) is programmed.

           else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
                       +(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
         end
       end
     end

The second and third end statements terminate the inner loops in i and j, so that all of the values of and are handled.

• Euler's method (Eq. (8)) is then used to advance the solution through nt=200 steps of length through the intermediate for loop in i2.

     w=w+wt*h;
     t=t+h;
   end

Note the use of the matrix facility of MATLAB (w and wt are arrays of dimensions nx = 21, ny = 21). The end completes the intermediate loop in i2.

• The numerical solution is displayed to produce the output in Table 5.

   fprintf('t = %5.2f   w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
 end

The end terminates the outer loop in i1.

• The analytical solution, Eq. (42), is then printed at the end (as reflected in Table 5).

 fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));

The programs in Listings 1, 2 and 3 have an essential feature that should be pointed out: only the operations of storing and retrieving numbers, arithmetic and looping are used (which are done naturally and very efficiently by computers). In other words, the integration of the three PDEs was not done by the usual analytical mathematics, but rather, was done essentially with arithmetic. This is the essence of the numerical method, that is, to replace mathematics of differentiation and integration with something much simpler, but performed naturally by computers with extraordinary speed and precision. This simplification was accomplished by replacing the partial derivatives in the PDEs with algebraic approximations, in this case, finite differences.

A central question then is how accurate the numerical solution is if approximations are used to compute it. We have observed the accuracy can be improved with more calculations (e.g., by reducing , and ). In other words there is a tradeoff between accuracy and calculations, and we basically assume the computer can do enough calculations to achieve the required accuracy. Historically, the increasing speed of computers has produced solutions to PDE problems with acceptable accuracy and continually increasing complexity. Thus, for the future, we can expect to compute numerical solutions of PDEs with essentially unlimited complexity if sufficient computing power is available.

References

External links