Boundary ConditionsOften, with PDEs, it is possible to determine a good numerical way to apply boundary conditions for a particular equation and boundary condition. The example given previously in the method of lines introduction section is such a case. However, the problem of finding a general algorithm is much more difficult and is complicated somewhat by the effect that boundary conditions can have on stiffness and overall stability. Periodic boundary conditions are particularly simple to deal with: periodic interpolation is used for the finite differences. Since pseudospectral approximations are accurate with uniform grids, solutions can often be found quite efficiently. Giving boundary conditions for partial differential equations. If you are solving for several functions then for any of the functions to have periodic boundary conditions, all of them must (the condition need only be specified for one function). If you are working with more than one spatial dimension, you can have periodic boundary conditions in some independent variable dimensions and not in others. This solves a generalization of the sineGordon equation to two spatial dimensions with periodic boundary conditions using a pseudospectral method. Without the pseudospectral method enabled by the periodicity, the problem could take much longer to solve. In[83]:= 
Out[83]=

In the InterpolatingFunction object returned as a solution, the ellipses in the notation {..., xmin, xmax, ...} is used to indicate that this dimension repeats periodically This makes a surface plot of a part of the solution derived from periodic continuation at t6. In[84]:= 
NDSolve uses two methods for nonperiodic boundary conditions. Both have their merits and drawbacks. The first method is to differentiate the boundary conditions with respect to the temporal variable and solve for the resulting differential equation(s) at the boundary. The second method is to discretize each boundary condition as it is. This typically results in an algebraic equation for the boundary solution component, so the equations must be solved with a DAE solver. This is controlled with the DifferentiateBoundaryConditions option to MethodOfLines. To see how the differentiation method works, consider again the simple example of the method of lines introduction section. In the first formulation, the Dirichlet boundary condition at x 0 was handled by differentiation with respect to t. The Neumann boundary condition was handled using the idea of reflection, which worked fine for a second order finite difference approximation, but does not generalize quite as easily to higher order (though it can be done easily for this problem by computing the entire reflection). The differentiation method, however, can be used for any order differences on the Neumann boundary condition at x 1. As an example, a solution to the problem will be developed using fourth order differences. This is a setting for the number of and spacing between spatial points. It is purposely set small so you can see the resulting equations. You can change it later to improve the accuracy of the approximations. In[85]:= 
This defines the vector of . In[86]:= 
Out[86]=

This discretizes the Neumann boundary condition at x 1 in the spatial direction. In[87]:= 
Out[87]=

This differentiates the discretized boundary condition with respect to t. In[88]:= 
Out[88]=

Technically, it is not necessary that the discretization of the boundary condition be done with the same difference order as the rest of the DE; in fact, since the error terms for the onesided derivatives are much larger, it may sometimes be desirable to increase the order near the boundaries. NDSolve does not do this because it is desirable that the difference order and the InterpolatingFunction interpolation order be consistent across the spatial direction. This is another way of generating the equations using NDSolve`FiniteDifferenceDerivative. The first and last will have to be replaced with the appropriate equations from the boundary conditions. In[89]:= 
Out[89]=

Now we can replace the first and last equation with the boundary condition. In[90]:= 
Out[92]=

NDSolve is capable of solving the system as is for the appropriate derivatives, so it is ready for the ODEs. In[93]:= 
Out[93]=

This shows a plot of how well the boundary condition at x 1 was satisfied In[94]:= 
Treating the boundary conditions as algebraic conditions saves a couple of steps in the processing. This replaces the first and last equations (from before) with algebraic conditions corresponding to the boundary conditions. In[95]:= 
Out[97]=

This solves the system of DAEs with NDSolve. In[98]:= 
Out[98]=

This shows how well the boundary condition was satisfied. In[99]:= 
For this example, the boundary condition was satisfied well within tolerances in both cases, but the differentiation method did very slightly better. This is not always true; in some cases, with the differentiation method, the boundary condition can experience cumulative drift since the error control in this case is only local. The Dirichlet boundary condition at x 0 in this example shows some drift. This makes a plot which compares how well the Dirichlet boundary condition at x 0 was satisfied with the two methods. The solution with the differentiated boundary condition is shown in black. In[100]:= 
When using NDSolve, it is easy to switch between the two methods by using the DifferentiateBoundaryConditions option. Remember that when you use DifferentiateBoundaryConditions>False, you are not as free to choose integration methods; the method needs to be a DAE solver. With systems of PDEs or equations with higher order derivatives having more complicated boundary conditions, both methods can be made to work in general. When there are multiple boundary conditions at one end, it may be necessary to attach some conditions to interior points. Here is an example of a PDE with two boundary conditions at each end of the spatial interval. This solves a differential equation with two boundary conditions at each end of the spatial interval. In[101]:= 
Out[101]=

Understanding the message about spatial error will be addressed in the next section. For now, ignore the message and consider the boundary conditions. This forms a list of InterpolatingFunctions differentiated to the same order as each of the boundary conditions. In[102]:= 
Out[102]=

This loads the package needed for logarithmic plots and makes a logarithmic plot of how well each of the four boundary conditions is satisfied by the solution computed with NDSolve. In[103]:= 
In[104]:= 
It is clear that the boundary conditions are satisfied to well within the tolerances allowed by AccuracyGoal and PrecisionGoal options. It is typical that conditions with higher order derivatives will not be satisfied as well as those with lower order derivatives. Inconsistent Boundary ConditionsIt is important that the boundary conditions you specify be consistent with both the initial condition and the PDE. If this is not the case, NDSolve will issue a message warning about the inconsistency. When this happens, the solution may not satisfy the boundary conditions, and in the worst cases, instability may appear. In this example for the heat equation, the boundary condition at x = 0 is clearly inconsistent with the initial condition. In[1]:= 
Out[1]=

This shows a plot of the solution at x = 0 as a function of t. The boundary condition is clearly not satisfied. In[2]:= 
Out[2]=

The reason the boundary condition is not satisfied is because once it is differentiated, it becomes , so the solution will be whatever constant value comes from the initial condition. When the boundary conditions are not differentiated, the DAE solver in effect modifies the initial conditions so that the boundary condition is satisfied. In[3]:= 
Out[4]=

It is not always the case that the DAE solver will find good initial conditions which lead to an effectively correct solution like this. A better way to handle this problem is to give an initial condition that is consistent with the boundary conditions, even if it is discontinuous. In this case the unit step function does what is needed: This uses a discontinuous initial condition to match the boundary condition, giving a solution correct to the resolution of the spatial discretization. In[5]:= 
Out[6]=

In general, with discontinuous initial conditions, spatial error estimates cannot be satisfied, since they are predicated on smoothness so, in general, it is best to choose how well you want to model the effect of the discontinuity by either giving a smooth function which approximates the discontinuity or by specifying explicitly the number of points to use in the spatial discretization. More detail on spatial error estimates and discretization is given in the next section. A more subtle inconsistency arises when the temporal variable has higher order derivatives and boundary conditions may be differentiated more than once. Consider the wave equation The initial condition satisfies the boundary conditions, so you might be surprised that NDSolve issues the NDSolve::ibcinc message. In this example, the boundary and initial conditions appear to be consistent at first glance, but actually have inconsistencies which show up under differentiation. In[8]:= 
Out[8]=

The inconsistency appears when you differentiate the second initial condition with respect to , giving , and differentiate the second boundary condition with respect to , giving = . These two are inconsistent at x = t = 0. Occasionally, NDSolve will issue the NDSolve::ibcinc message warning about inconsistent boundary conditions when they are actually consistent. This happens due to discretization error in approximating Neumann boundary conditions or any boundary condition which involves a spatial derivative. The reason this happens is that spatial error estimates (see the next section) used to determine how many points to discretize with are based on the PDE and the initial condition, but not the boundary conditions. The onesided finite difference formulas which are used to approximate the boundary conditions also have larger error than a centered formula of the same order, leading to additional discretization error at the boundary. Typically this is not a problem, but it is possible to construct examples where it does occur. In this example, because of discretization error, NDSolve incorrectly warns about inconsistent boundary conditions. In[9]:= 
Out[9]=

A plot of the boundary condition shows that the error, while not large, is outside of the default tolerances. In[10]:= 
Out[10]=

When the boundary conditions are consistent, a way to correct this error is to specify that NDSolve use a finer spatial discretization. With a finer spatial discretization, there is no message and the boundary condition is satisfied better. In[11]:= 
Out[12]=

