Introduction

Most physical phenomena in fluid dynamics, electricity, electromagnetism, mechanics, classical optics or in heat flow are described by partial differential equations (PDEs). In fact, well-known laws of physics, such as Maxwell’s equations, the Navier–Stokes equations, the heat equation, the wave equation and Schrödinger’s equation of quantum mechanics, are stated in terms of PDEs; that is, these laws describe physical phenomena by relating space and time derivatives. Derivatives in these equations represent quantities like velocity, acceleration, force, friction, flux and current. Hence, we have equations relating partial derivatives of some unknown quantity that we would like to find [1–3].

A PDE is an equation containing various partial derivatives of a multivariable function. Unlike ordinary differential equations (ODEs), where the unknown function depends only on one variable, the unknown function depends on several variables in PDEs. For instance, a temperature field may depend on both location x and time t [2].

Throughout this monograph, the following notational convention of partial derivatives will be used for simplicity:

With this abbreviated notation at hand, we list a few well-known PDEs [3] with unit coefficients:

The unknown function u is called the dependent variable and is a function of several independent variables like the time t and spatial coordinates like x. Partial differential equations are constructed from terms where the dependent variables are differentiated with respect to the independent variables. For example, in the one-dimensional heat equation , the dependent variable is a function of two independent variables x and t, and in the Laplace equation in polar coordinates , the dependent variable depends on r and θ [4].

A function that satisfies a PDE is called a solution of the PDE [5] and is typically what is unknown and sought.

Since the general theory and solution methods usually apply only to a given class of equations, a classification of PDEs is important. Basic classifications of PDEs can be listed as [6]:

1.   Order of the PDE . The order of a PDE is the order of the highest partial derivative in the equation, for example:

2.   Number of Variables . The number of variables is the number of independent variables, for example:

The most general first-order PDE in two independent variables can be written as:

where denotes a function of its arguments, i.e. the highest order of derivatives in is one.

Similarly, the most general second-order PDE in two independent variables can be written as:

3.   Kinds of Coefficients . If the coefficients , , , , and in equation ( 7 ) are constants, then ( 8 ) is called a PDE with constant coefficients; otherwise, it is called a PDE with variable coefficients [ 9 ].

4.  Linearity. Partial differential equations could be either linear or nonlinear. If the dependent variable u and all its partial derivatives occur linearly in the PDE, then the PDE is linear. More precisely, a second-order linear PDE in two independent variables is an equation of the form

where , , , , , and can be constants or given functions of independent variables only and not all , , are zero simultaneously [10].

Due to the importance of the concept, let us give a few more examples:

In general, if is a differential operator satisfying and for any functions , and any constant , then is a linear differential equation .

Consider the given PDE ; we can see that is linear since:

Similarly, for the PDE , we can see that is nonlinear since:

An important property that results from linearity of the differential operator is that if and are both solutions, so is . This is also called the superposition principle .

5.   Homogeneity . Equation ( 11 ) is called homogeneous if the right-hand side is identically zero. Otherwise, i.e. if is not identically zero, then the equation is called inhomogeneous [ 12 ]. Another consequence of linearity is that with being a linear differential operator, if a solution of (homogeneous solution) is added to a solution of (inhomogeneous solution), then an inhomogeneous solution is obtained [ 13 ].

In this monograph, we will be studying PDEs following the previous classification. We will start with the simplest possible PDE and then continue with more difficult ones. For each class, we will be going over the methods that are available to solve those problems. Note that for PDEs for which symbolic solution methods are yet not available, we recommend using NDSolve.

Recall that the general solution of a linear homogeneous ordinary differential equation (ODE) of order m is a linear combination of m independent solutions with m arbitrary constants [14]. Let us see what changes in this regard when we are dealing with PDEs in simple cases.

The simplest possible PDE is (where ). Its general solution is , where is any function of one variable. For instance, , , and are several solutions. Since the solutions do not depend on , they are constant on the lines in the xy plane, where is a constant.

For this PDE, DSolve returns:

The following Manipulate shows that the solutions are constant on the curves:

Suppose we want to find solutions of . If we integrate the PDE with respect to x once, we get , where is an arbitrary function. Since the PDE involved a variable y in addition to the variable x, the integration constant became instead of a regular constant coefficient. Integrating with respect to x again leads us to . Note that there are two arbitrary functions in the solution, which play the same role as integration constants do in ODEs [15].

DSolve returns an equivalent result for this PDE:

Suppose we want to find solutions of . In fact, since there is no differentiation with respect to y and there is no explicit appearance of y in the PDE, this equation is really an ODE with an extra independent variable y. Since a general solution to the ODE is , using arbitrary functions and as integration constants instead of and , the solution to the PDE is [16].

For this PDE, DSolve returns the same result:

Suppose we want to find solutions of . First, we integrate the PDE with respect to x while treating y as fixed and we get , where is an arbitrary function. Next, we integrate with respect to y while treating x as fixed, and we get the solution , where and (such that ) are arbitrary functions [17].

DSolve returns the same result for this PDE:

Initial and Boundary Conditions

As we saw in previous examples, PDEs typically have solutions involving arbitrary functions. In order to single out a specific solution, additional auxiliary conditions are imposed. In other words, we impose the conditions to specify a unique solution. These conditions are motivated by the nature of the problem and they come as initial conditions and boundary conditions [18].

An initial condition defines the physical state of a time-dependent PDE at a precise time .

For instance, for the diffusion equation, an initial condition is

where is a known function with vector denoting the space variables. For a diffusing substance, is the initial concentration, whereas for heat flow, is the initial temperature [19].

For the wave equation, initial conditions are given in pairs:

where is the initial position and is the initial velocity. It is essential from the physics of the problem that both of them have to be specified in order to determine the position uniquely at later times [20].

In summary, for each order of time derivative of the PDE, one initial condition has to be specified.

PDEs that come from real-world applications are limited to be operational in a domain Ω. For instance, for a vibrating string with length l, Ω is the interval . The boundary of Ω, denoted by ∂Ω, involves only the two points and . In the case of a drumhead, the domain is a plane region and its boundary is a closed curve surrounding the region, whereas for the diffusing chemical substance, Ω is the container holding the liquid, and hence ∂Ω is its boundary surface [21].

Boundary conditions are the interface between the partial differential equation operational in the domain Ω and everything that is outside of that domain and not part of the solution. It is necessary to specify some boundary conditions for a unique solution to be determined. The three most important and common types of boundary conditions are:

1.  Dirichlet condition: is specified on the boundary ∂Ω.

2.  Neumann condition: the normal derivative is specified on the boundary ∂Ω.

3.  Robin condition: is specified on the boundary ∂Ω.

In these conditions, a can be a function of x and t, and each condition is supposed to hold for all t and for all x on the boundary [22].

Note that although there is a relationship between the system function NeumannValue and the Neumann or Robin conditions as used here, their definitions are different.

In Figure 1, locations where Dirichlet conditions might be specified are shown in blue on the boundary of the domain Ω. Dirichlet conditions specify the solution values satisfied at those points. Locations where Neumann conditions might be specified are shown in green on the boundary of the region Ω and they specify the inward flux in the direction of the outward normal on the boundary .

Usually, we write boundary conditions as equations. For instance, a Neumann condition is written as the equation

where g is a given function. Any of these boundary conditions are called homogeneous if the specified function . Otherwise, they are called inhomogeneous. In Figure 1, denotes the unit normal vector on ∂Ω, pointing outward from Ω, and denotes the directional derivative of u in the outward normal direction [23].

In one-dimensional problems where Ω is just an interval , the boundary consists of just the two endpoints, and these boundary conditions take the simple form

where , , and can be functions of t. It is possible to have different types of boundary conditions at different locations on the boundary [24].

A time-dependent PDE that is specified with an initial condition is called an initial value problem (IVP), whereas a PDE with boundary conditions is called a boundary value problem (BVP). A time-dependent PDE that involves both initial and boundary conditions is called an initial value-boundary value problem (IVBVP).

Suppose we want to solve the differential equation with a boundary condition given as , where and .

We have seen in an earlier example that the equation has the general solution . If we now impose the given boundary condition , we find that , and hence the solution to the problem is uniquely determined as .

For this BVP, DSolve returns:

Suppose we want to solve the following problem with constant coefficients [25]:

where and .

For this problem, DSolveValue returns:

While solving the preceding problem, DSolveValue uses a method that is usually referred to as the "Fokas method" or the "uniform transform method". Unlike the classical methods for solving linear partial differential equations with constant coefficients that rely on separation of variables and specific integral transforms, the Fokas transformation method is not limited to specific equations with special boundary conditions. The details of the method are beyond the scope of the monograph, and hence we refer to [26] for details.

First-Order Linear Equations—Method of Characteristics

The general first-order linear PDE IVP with two independent variables is given as:

One solution technique to solve first-order linear PDEs is the method of characteristics, where we aim to find a change of independent variables to new variables in order to obtain an ODE IVP that is easier to solve than (27) [28].

Before getting into the details, to motivate the method, suppose that we know the solution to problem (29). Suppose also that the solution is smooth and we have a smooth parametric curve, which we will call a characteristic curve, in the x-t plane such that and . This means that and are functions of a single parameter . Now consider how our solution would change as we move along these parametric curves:

Notice that we can use the similarity between the preceding relationship (30) and the left-hand side of the PDE in (31). If we can find and such that:

the PDE in (32) will reduce to the ODE:

Hence, with the knowledge of , i.e. the initial point of the parametric curve, is determined uniquely.

At this point we have to answer the question: How do we find the characteristic curves that reduce the problem (33) into an ODE IVP?

In order to find the characteristic curve and solve the resulting ODE (34) uniquely, we need to define an IVP. We require that new coordinates s (to parameterize x and t) and τ (for the initial curve) have the following properties [35]:

1.   varies along the characteristic curves , ranging from 0 to starting at the initial curve, the line for this case.

2.   varies along the initial curve and stays constant along each characteristic curve.

The following figure illustrates these properties [36]:

Hence to uniquely define the characteristic curves in the new variables s and τ, we define the IVPs:

and

Since , we have when . Hence in the new coordinate system, is mapped onto . Therefore, the initial condition is transformed into to act as an initial condition for the reduced ODE (37) in s.

By applying the chain rule and the preceding properties, we arrive at the ODE IVP:

So along the characteristics curves , the PDE problem (38) is reduced down to an ODE IVP (39). Once we solve this newly constructed problem, we change back to the original independent variables to obtain the solution to the original PDE (40) [41].

In summary, this approach, where:

3.  New variables define characteristic curves , each of which has as the initial point, and the initial data propagates along these curves,

4.   The s -coordinate indicates points along a characteristic curve,

5.   The -coordinate indicates the initial point ( ) for a given curve,

6.   For each , the initial value is evolved along the curve starting at the point according to the ODE obtained by expressing the PDE in characteristic variables,

is known as the method of characteristics.

The characteristic equations for (42), , , are autonomous. The term autonomous is used when the equation has no explicit dependence on s.

The characteristics satisfy the ODE:

Hence for problem (43), the characteristics can be solved independently of the solution u.

Equation (44) is a fundamental equation used when applying this method, as will be shown in the following examples.

Note that if in (45), then u remains constant along the characteristic curves.

Let us solve the constant coefficient equation

where a and b are constants and not both zero. Having y take the role of t and using (46), the characteristic curves satisfy the ODE . Solving this ODE gives:

are the characteristics. Again, remains constant on each such curve. So , where f is an arbitrary function. Therefore, the general solution of (47) is obtained by solving for . That is:

Hence

and we get an equivalent solution directly by DSolve:

We can easily verify the preceding solution:

The equation

is linear and homogeneous but has a variable coefficient (y), for which the characteristic curves satisfy the ODE . Solving the ODE gives:

are the characteristics. Again, is a constant on each such curve. So , where f is an arbitrary function. Therefore, the general solution of (48) is obtained by solving for . That is:

Hence

We get the same solution directly from DSolve:

Geometrically, the "picture" of the solution is that it is constant on each characteristic curve:

Find the solution of (49) that satisfies the auxiliary condition . Indeed, putting in (50), we get , so that . Therefore, . The same solution can be found by DSolve with the call:

Solve the PDE

The characteristic curves satisfy the ODE . Solving the ODE gives:

are the characteristics. Again, is a constant on each such curve. So , where f is an arbitrary function. Therefore, the general solution of (51) is obtained by solving for . That is:

Hence

and we get the same solution directly by DSolve:

Suppose we want to solve the IVP with constant coefficients [52]:

First, to find the characteristic curves, solve the following ODE IVPs:

for which DSolve gives:

Hence, our characteristic curves are:

Using this new coordinate system defined by (53), we reduce the PDE to the ODE:

Solving this IVP, we get:

So

is the solution in the sτ-coordinate system. In order to find u as a function of x and t, solve the transformation

for s and τ in terms of x and t. We have:

Hence, our answer is

which we could have directly found in an equivalent form by DSolveValue:

If we plot the solution at three different time values, we get:

The preceding plot shows the initial wave in orange. Note that with time progressing to and , the wave dampens out; in other words, the amplitude decreases. Also, note that the peaks at the various times move slightly to the right.

Nonlinear First-Order Equations—Conservation Equations

As an initial example, we consider the problem of modeling traffic flow. Let represent the density of cars at the point x and at the time t of a highway with no exit or entrance ramps. Then on any segment of the highway, u satisfies the following equation [54]:

where is the flux, i.e. the cars per minute passing the point x. Both sides of the preceding equation represent the change in the number of cars within .

Using the fundamental theorem of calculus, we obtain:

which is rearranged to:

Because the segment is arbitrary, we conclude that u is a solution of the one-dimensional conservation equation:

We now proceed to solve the conservation equation IVP:

Using the chain rule, we rewrite the PDE as:

or

where .

So the IVP that we would like to solve is:

Similar to the linear case, making use of the method of characteristics, the characteristics satisfy the equations:

The PDE reduces to the ODE:

with initial condition

DSolve tells us that the characteristic curves are defined by:

and:

The ODE IVP for u has the solution:

Hence, since

it follows that

and the original IVP (55) has the implicitly defined solution:

Suppose we want to solve the following IVP:

Since and , the implicit solution of the IVP is:

which is the solution found by DSolveValue:

In this particular example, although we cannot do this in many cases, we can solve for explicitly in terms of x and t.

Suppose we want to solve the PDE:

Rewriting the equation as:

with , the implicit solution of the PDE is:

which is the solution found by DSolve:

Here is the verification of the implicit solution our PDE:

Well-Posed Problems

A well-posed problem is a PDE and a set of initial and/or boundary conditions (or other auxiliary conditions) in a domain with the following fundamental properties [56]:

1.   Existence , that there exists at least one solution satisfying all the PDE and the conditions, and

2.   Uniqueness , that there is at most one solution, and

3.   Stability , that the unique solution varies continuously on the data of the problem; i.e. if the data is changed a little, the solution changes a little as well.

When a physical problem is modeled by a PDE, it is expected that we formulate the problem with physically realistic auxiliary conditions that make the problem well-posed. In the case that too few auxiliary conditions are imposed, there may be more than one solution (non-uniqueness) and the problem is called underdetermined. However, if there are too many auxiliary conditions, there may be no solution at all (non-existence) and the problem is called overdetermined.

To illustrate the concept of a well-posed PDE, consider the following setup. A vibrating string with an external force, whose ends are moved in a specified way, satisfies the problem [57]:

for . The five functions , , , and form the data for this problem. Existence and uniqueness mean that there is exactly one solution for differentiable functions f, φ, ψ, g and h. Stability means that if any of these five functions are slightly changed, then u also gets changed only slightly.

Types of Second-Order Linear Equations

A second-order linear equation in two independent variables is an equation of the form [58, 59]:

where , , , , , and can be constants or given functions of x and y.

The independent variables and in (60) are usually used to denote space variables. In the case that we have the time variable instead of , then we will have terms , and in (61). From the classification point of view, this does not make any difference; only physical meaning changes.

There are three basic types of second-order linear equations of the form (62):

1.   Parabolic equations (if , reducible to the canonical form , describe heat flow and diffusion processes and satisfy the property),

2.   Hyperbolic equations (if , reducible to the canonical forms or , describe vibrating systems and wave motion),

3.   Elliptic equations (if , reducible to the canonical form , describe steady-state phenomena and satisfy the property),

where denotes terms of order 1 or 0. This will be clarified by the following examples.

The Heat Equation

In this section, we consider heat flow problems starting with a simple model and build further on this. We refer to the tutorial on heat transfer models for a general introduction.

Fixed End Temperatures and Separation of Variables

Suppose that we have a rod with length l whose lateral sides are wrapped with insulation, so that heat can flow in and out of the rod at the ends, but not across the lateral boundary. If this rod is placed in an environment whose temperature is fixed at some temperature , after a sufficient time, the temperature of the entire rod stabilizes and comes to a steady-state temperature similar to the environment. Then suppose we take the rod out of the environment at a time that we call and attach two temperature components to the ends of the rod to keep the ends at specific fixed temperatures, say and . If we monitor the temperature chart of the rod, the basic PDE that is modeling the one-dimensional heat flow is:

PDE (63) relates the quantities and , where the former represents the rate of change in temperature with respect to time, and the latter represents the concavity of the temperature profile . In other words, PDE (64) compares the temperature at one point to the temperature at neighboring points. The proportionality constant k is a property of the material that the rod is made of [65, 66].

Since the temperature was fixed for all times at and at the two endpoints and , for the boundary conditions, we could say:

Recall that we started monitoring the rod temperature from the time the rod had reached a temperature of ; we have the initial condition:

Hence the PDE (67), Dirichlet BCs (68) and the IC (69) make up the IVBVP.

Let us first suppose that , i.e. the ends of the rod are held at zero temperature. So, we want to solve the heat flow IVBVP:

We assume that the solution of (70) can be written as:

i.e. the solution is separable.

Substituting (71) into (72) leads to:

Dividing both sides of (73) by results in:

where has to be constant since is a function of solely x and is a function of solely t.

Equation (74) leads to two ODEs:

for which the general solutions are:

where and are constants.

The next step is imposing (75) on (76), which requires that . Thus we have:

Since we are interested in a nontrivial solution, i.e. and should not both be zero, we arrive at:

So for ,

are distinct solutions of (77).

Hence there are infinitely many solutions of (78) such that:

where are constants. Due to the superposition principle for linear equations, the sum of solutions is also a solution. Therefore,

is also a solution of (79) and (80).

Imposing (81) on (82), we get:

In the preceding, are called eigenvalues and are called eigenfunctions.

Using the orthogonality of the eigenfunctions, i.e.:

and:

we find that:

Substituting (83) into (84), we conclude that:

which is the same solution that DSolveValue returns:

The method outlined is the so-called separation of variables method, also called the eigenfunction expansion method.

Note that the and cases do not lead to nontrivial solutions, and hence we omitted these cases.

Now let us return to the IVBVP formed by the PDE (85), Dirichlet BCs (86) and the IC (87). For this IVBVP, DSolveValue returns:

For , , , and , a three-dimensional plot of the solution is:

Notice that while plotting the solution we truncated the infinite summation at to get a sufficiently good approximation for the summation by using the function TruncateSum. Clearly, if we truncate at larger values of , we could get better approximations. In this example, truncating at was sufficient for our purpose.

In the next subsection, we will see that the first two terms in the preceding solution form the steady-state solution.

Steady-State Temperatures

Suppose that the function represents temperature. After a sufficiently long time under the same environment and conditions, experiences show that the variation of temperature with time dies away. Hence, we expect that the limit of , as t tends to infinity, exists and depends only on x [88]:

and also that

The steady-state temperature distribution function has to satisfy the BCs (89) and the heat equation (90), which are valid for all . Therefore, should be the solution to the problem:

which is equal to the limit of the solution of the preceding IVBVP as .

Let , , , and in the preceding IVBVP and solve the problem with concrete values:

The following Manipulate computation shows that after a long enough time, the solution starts converging to the steady-state solution:

Note that Quiet in the computation is used to suppress General::munfl messages that get generated during the computation.

Insulated Rod

Suppose now that the ends of the rod at and are insulated instead of being held at constant temperatures, i.e. Neumann BCs are given. The IVBVP that describes the temperature in this rod is [91]:

For this IVBVP, DSolveValue returns the solution:

Let , and in the preceding IVBVP, and solve the problem with concrete values:

If we truncate the infinite summation at, for example, , the three-dimensional plot of the solution is:

Nonhomogeneous Equation with Heat Source

The nonhomogeneous equation

models the case where the rod has an internal heat source everywhere along the rod and for all time t [92]. With the same set of BCs (93) and the IC (94) from an earlier example, in this case DSolveValue returns:

Let , , , , , in the preceding IVBVP, and solve the problem with concrete values:

If we truncate the infinite summation at, for example, , the three-dimensional plot of the solution is:

In the case that a wire carrying electrical current passes through the rod and the resistance generates a constant heat source , DSolveValue then returns [95]:

Let , , , , , in the preceding IVBVP, and solve the problem with concrete values:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution for the problem with constant heat source is:

Boundary Conditions for Diffusion-Type Problems

In this subsection, we consider examples of IVBVPs with conditions on the boundary given as Dirichlet conditions, Neumann conditions or Robin conditions.

Suppose that for the heat-flow problem over a one-dimensional rod (96) with , the temperatures on the boundary are specified, i.e. the conditions at both ends are given as the Dirichlet conditions, and IC is given as [97]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Note that IC, :

and the BCs, and :

are satisfied by the solution.

Suppose that for the heat-flow problem with a source over a one-dimensional rod (98), where and , the fluxes of the temperature at the boundary are specified, i.e. the conditions at both ends are given as the Neumann conditions, and the IC is given as [99]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that for the heat-flow problem with a source over a one-dimensional rod (100), where and , the temperature at one end of the boundary and the flux of the temperature at the other end are specified, i.e. the condition at one end, namely at , is given as the Dirichlet condition, and at the other end, namely at , is given as the Neumann condition. Suppose also that the IC is given as ; then we have [101]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that for the heat-flow problem with a source over a one-dimensional rod (102), where and , the temperature at one end of the boundary and the temperature of the surrounding medium at the other end are specified, i.e. the condition at one end, namely at , is given as the Dirichlet condition, and at the other end, namely at , is given as the Robin condition. In this problem at , the end of the bar is in contact with a surrounding medium with a prescribed temperature [103]. In other words, the right end is exposed to convective heat transfer [104]. Suppose also that the IC is given as ; then we have:

The preceding solution is a result of the method of separation of variables and hence the infinite summation involves summation over eigenfunctions corresponding to eigenvalues . The eigenvalues are defined by the condition [105]:

In order to plot the solution, we compute the eigenvalues subject to to get a sufficiently good approximation for the infinite summation in the solution:

So there are 503 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 503 eigenvalues:

A three-dimensional plot of the truncated solution is:

Suppose that for the heat-flow problem with a source over a one-dimensional rod (106), where and , the flux of the temperature at one end is specified and the other end is in contact with another medium, i.e. the condition at one end, namely at , is given as the Neumann condition, and at the other end, namely at , is given as the Robin condition. Suppose also that the IC is given as ; then we have [107]:

Let us compute the eigenvalues s.t. :

So there are 504 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 504 eigenvalues:

A three-dimensional plot of the truncated solution is:

Suppose that for the heat-flow problem with a source over a one-dimensional rod (108), where and , both ends are in contact with another medium [109], i.e. the BCs at both ends are given as Robin conditions. Suppose also that the IC is given as ; then we have:

Let us compute the eigenvalues such that :

So there are 504 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 504 eigenvalues:

A three-dimensional plot of the truncated solution is:

Heat Loss into the Surrounding Medium

Suppose that a uniform rod is surrounded by a medium with fixed temperatures at both ends, and it encounters heat loss through its lateral surface into its surrounding medium. Then the heat flow PDE describing this situation is:

where the term is for the heat loss across the lateral boundary [110].

For , and , if both the left and right ends of the rod are held at the fixed temperature zero, and if the initial temperature distribution is given as , then we have the following heat-flow problem:

For this IVBVP, DSolveValue returns:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Without the heat loss term with the same initial and boundary conditions, DSolveValue returns:

To see the effect of the heat loss term, let us plot the solutions side by side:

Diffusion Problem in the Semi-infinite Rod

In all of the preceding problems, we have worked over finite intervals. If the rod under consideration is very long, we can treat it as semi-infinite, i.e. , and if the rod is uniform and there is no external source generating heat, the partial differential equation describing the temperature remains [111, 112]:

Suppose that at the temperature is held constant, i.e. . Due to the absence of another boundary, there is no other boundary condition. However, we expect that remains finite, i.e. bounded, as . Hence our model is:

For , and , DSolveValue returns:

Assumptions is an option for DSolveValue, and in the preceding computation it is used to specify the domain for the independent variable x.

A three-dimensional plot of the solution is:

Suppose that the BC at is given as a Neumann condition:

A three-dimensional plot of the solution is:

Diffusion Problem in the Infinite Rod

Suppose that we want to study heat conduction in the center of a very long rod. In this case, we can assume that it extends from to . Hence there are no boundary conditions, and the problem to be solved is [113]:

For and , DSolveValue returns:

A three-dimensional plot of the solution is:

Suppose there is a source that is generating heat in the system:

A three-dimensional plot of the solution is:

Transforming Hard Equations into Easier Ones

The diffusion-convection equation [114]:

can be transformed into the well-studied heat equation

by using the transformation

Of course, the IC and the BCs in the problem must be transformed too using the same transformation. Suppose that we have an IVBVP with , where the IC is given s.t. and Dirichlet BCs at both ends are given s.t. . For this IVBVP, DSolveValue returns:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Fourier Transform and the Heat Equation

Consider the heat flow in an infinite rod where the initial temperature is known, and suppose that we want to solve the IVP:

which is also called a Cauchy problem [115].

First, let , and transform the equation:

Then transform the IC, :

Hence we have an ODE in with an initial condition. Solving this ODE-IVP, we get:

Now, transforming the preceding solution back gives the solution to the original PDE-IVP:

Applying ExpToTrig to this solution and simplifying gives:

The same solution could be found directly by DSolveValue:

A three-dimensional plot of the solution is:

Laplace Transform and the Heat Equation

Imagine a cylindrical rod of length with perfectly insulated curved sides so that no heat can enter or escape. At time , the temperature of the rod is given by , . For , temperatures at both ends of the rod are given as and . The IVBVP for temperature in the rod is [116]:

First, let , and transform the PDE into an ODE using LaplaceTransform:

Note that the IC, , is used in the preceding ReplaceAll.

Transforming the BCs gives:

Solving the ODE-IVP gives:

Now, transforming the preceding solution back gives the solution to the original IVBVP:

The same solution could be found directly by DSolveValue:

A three-dimensional plot of the solution is:

Two-Dimensional Heat Equation

Suppose that we study the temperature distribution in a thin rectangular plate over the finite two-dimensional domain where and , for which the lateral surfaces of the plate are insulated, so there is no heat loss through the lateral surfaces. The boundaries and are kept at a fixed zero temperature (Dirichlet conditions). The boundary is insulated (Neumann condition), and the boundary is in contact with a surrounding medium with temperature zero (Robin condition). The system has no sources that generate heat. The initial temperature distribution is given as , and the constant thermal characteristics of the plate material are . Hence the IVBVP that we would like to solve is [117]:

For this IVBVP, DSolveValue finds the solution:

Let us compute the eigenvalues such that :

So there are 100 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 100 eigenvalues:

Now we visualize the three-dimensional plot of the solution:

Three-Dimensional Heat Equation

Suppose that we study the temperature distribution in a solid in the form of a rectangular parallelepiped where , and . The boundaries , , and are kept at a fixed temperature of zero degrees. The boundary is insulated, and the boundary is in contact with a surrounding medium at temperature zero. The system has no sources that generate heat. The initial temperature distribution is given as , and the constant thermal characteristics of the solid material are . Hence the IVBVP that we would like to solve is [118]:

For this IVBVP, DSolveValue finds the solution:

Let us compute the eigenvalues such that :

So there are 202 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 202 eigenvalues and truncate the infinite sums at, for example, :

A three-dimensional contour plot of the solution for :

Heat Flow IVBVPs in Polar and Cylindrical Coordinates

In this subsection, we study heat flow IVBVPs in circular regions and cylindrical domains. In the circular region, we first consider the case with circular symmetry, i.e. no angular dependency. In the second example, we consider the case with angular dependency. The third example deals with the heat-flow equation in a cylindrical domain.

Suppose that we study the heat conduction in a thin circular plate with radius r. Suppose also that the system is circularly symmetric, thermal coefficients are spatially invariant and there are no heat-generating sources. Then the PDE that we are dealing with is [119]:

Suppose we want the solutions to this partial differential equation over the finite interval ; then we require boundedness for , i.e. ∃_M such that ≤M, ∀_tt>0. With , suppose that we have the following BCs:

and the following IC:

DSolveValue returns the following solution to this IVBVP:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Suppose that we study the temperature distribution of a thin circular plate with radius r over the two-dimensional region and , with no heat-generating sources. Suppose also that the circular plate is very thin. Hence z-dependence is neglected and the PDE that we are dealing with is [120]:

One of the BCs is that we require boundedness for . With , suppose we have the following BCs:

and the following IC:

DSolveValue returns the following solution to this IVBVP:

For practical purposes, we should set in the preceding solution. Hence we need to compute the first 50 eigenvalues for for :

Now compute the summation using these 50 eigenvalues:

We visualize the three-dimensional plot of the solution:

Now we visualize the three-dimensional plot of the solution in the cylindrical domain:

Suppose that we study the temperature distribution of a cylinder with radius r over the three-dimensional region , and . Hence the PDE that we are dealing with is [121]:

One of the BCs is that we require boundedness for . With , suppose that we have the following BCs:

and the following IC:

DSolveValue returns the following solution to this IVBVP:

For practical purposes, we should truncate the infinite summation at, for example, :

We visualize the three-dimensional contour plot of the solution for :

The Wave Equation

In this section we consider wave phenomena problems starting with a simple problem and build further on this. We refer to the tutorial on wave equation for a general introduction.

One-Dimensional Wave Equation

Suppose that we consider the small vibrations of a string with length l. We assume the string is fixed at each end, and the vibrations of the string occur in a plane, i.e. the displacement is parallel to the string. If we monitor the longitudinal displacement of the string from equilibrium, denoted by , the basic PDE that is modeling the longitudinal and torsional vibrations of a string is the wave equation [122]:

where c is a constant physical parameter related to the linear density of the material [123], also called the wave speed [124]. The term represents the net force due to the tension of the string. The wave equation contains a second-order time derivative , and hence to define the solution uniquely for , it requires two ICs:

which is different from the heat equation problem case, where only one IC was required [125].

As an example, suppose that we have a string with length 1 and ; then the following PDE, BCs and ICs make up the IVBVP:

for which DSolveValue returns:

and if we truncate the infinite summation at, for example, , a two-dimensional plot of the solution for various values of x is as follows:

D'Alembert's Solution

It is possible to articulate the solution of the wave equation (126) directly in terms of the initial data in some cases. Let , and . Then the wave equation (127) becomes [128]:

and DSolve returns the following solution for this new PDE:

Hence if we now transform back to the original variables, the solution to the wave equation (129) is:

An equivalent form of the solution can be found by DSolveValue directly:

Next, if we impose initial conditions on the wave equation on the infinitely long string and seek the solution, we get the solution in terms of the initial data:

This solution is known as d'Alembert's solution or the traveling wave solution, representing the solution as the superposition of two waves, one moving to the left and the other to the right, with propagation speed c. This form of the solution shows how the initial data affects the solution at later times [130].

Generalization of the One-Dimensional Wave Equation: The Telegraph Equation

Suppose this time that we have a damping effect in our model as well as an external force acting on the system. Suppose also that we have a restoring force that is directed opposite to the displacement of the string. Hence the PDE modeling the wave phenomena becomes:

where denotes the spatial-time-dependent wave amplitude, c denotes the wave speed, β denotes the restoring force coefficient, γ denotes the damping coefficient of the medium and denotes the external force. The PDE (131) is also known as the telegraph equation. Assuming that the medium is uniform, all the preceding coefficients are spatially invariant and can be written as constants [132, 133].

Assume that

then DSolveValue returns:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Inhomogeneous Boundary Conditions

Let us consider a problem with sources given at the boundary, i.e. inhomogeneous BCs are given. So suppose that we have an IVBVP such as [134]:

For this IVBVP, DSolveValue returns:

Truncating the infinite summation at, for example, , we get:

A three-dimensional plot of the truncated solution is:

Wave Equation in Unbounded Regions

In this subsection, we consider wave phenomena problems over the half-line and over the whole line.

Consider the problem:

over the half-line , and DSolveValue returns [135]:

A three-dimensional plot of the solution follows as:

Consider the problem:

over the whole line , and DSolveValue returns [136]:

A three-dimensional plot of the solution follows as:

Consider the nonhomogeneous IVP:

over the whole line , and DSolveValue returns [137]:

Let us make a list of data points from the solution to use in ListPlot:

Using ListPlot, we get the following three-dimensional plot:

Boundary Conditions Associated with the Wave Equation

In this subsection, we consider wave phenomena IVBVPs with Dirichlet, Neumann and Robin BCs.

Suppose that for the vibration problem over a finite string, we have controlled endpoints, i.e. the conditions at both ends are given as the Dirichlet conditions [138]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that for the vibration problem over a finite string, forces are specified on the boundaries such that the ends of the string are allowed to slide vertically on frictionless sleeves, i.e. the conditions at both ends are given as the Neumann conditions [139]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that for the vibration problem over a finite string, we have a controlled end and a force is specified on the other end, i.e. the condition at one end is given as the Dirichlet condition and at the other end is given as the Neumann condition [140]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that for the vibration problem over a finite string, we have a controlled end and an elastic attachment on the other end, i.e. the condition at one end is given as the Dirichlet condition, and at the other end is given as the Robin condition [141]:

Let us compute the eigenvalues such that :

So there are 202 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 202 eigenvalues:

A three-dimensional plot of the truncated solution is:

Suppose that for the vibration problem over a finite string, we have elastic attachments on both ends, i.e. the conditions at both ends are given as the Robin conditions [142]:

Let us compute the eigenvalues such that :

So there are 403 eigenvalues such that .

Next, let us compute the summation in the preceding solution corresponding to the first 403 eigenvalues:

A three-dimensional plot of the truncated solution is:

Fourier Transform and the Wave Equation

Suppose we want to solve the following IVP describing a vibration of an infinitely long string [143]:

First, let and transform the equation:

Then transform the ICs:

Hence we have an ODE in with an initial condition. Solving this ODE-IVP, we get:

Now, transforming the preceding solution back gives the solution to the original PDE-IVP:

Applying ExpToTrig to this solution and simplifying gives:

The same solution could be found directly by DSolveValue:

A three-dimensional plot of the solution is:

Laplace Transform and the Wave Equation

Suppose we want to solve the following IVBVP describing a vibration [144]:

First, let , and transform the PDE into an ODE using LaplaceTransform:

Note that the ICs, and , are used in the preceding ReplaceAll.

Transforming the BCs gives:

Solving the ODE-IVP gives:

Now, transforming the previous solution back gives the solution to the original IVBVP:

The same solution could be found directly by DSolveValue:

A three-dimensional plot of the solution is:

Two-Dimensional Wave Equation

Suppose that we study the transversal vibrations in a thin rectangular membrane over the finite two-dimensional rectangular domain where and . We assume that the membrane is in a medium with no damping. We also assume that the boundaries are held fixed, and the membrane has an initial displacement of and an initial speed of . There are no external forces acting on the system, and the wave speed is . Hence the IVBVP that we would like to solve representing the amplitude of the transversal wave is [145, 146]:

For this IVBVP, DSolveValue finds the solution:

Let us truncate the infinite summations at, for example, :

Now we visualize the three-dimensional plot of the solution as:

Three-Dimensional Wave Equation

Suppose that we study the vibrations in a rectangular parallelepiped where , and . We assume that there is no damping in the medium. We also assume that the boundaries are held fixed, and the object has an initial displacement of and an initial speed of . There are no external forces acting on the system, and the wave speed is . We seek the wave amplitude . Hence the IVBVP that we would like to solve is [147, 148]:

For this IVBVP, DSolveValue finds the solution:

Let us truncate the infinite summations at, for example, .

For :

A three-dimensional contour plot of the solution is:

Symmetric Vibrations of a Circular Membrane

Suppose that we consider the problem describing the displacement of a circular membrane with radius a that is fixed at its edge. Let us first treat the simple case in which the initial conditions are independent of θ. The wave equation [149]:

defining the displacement becomes:

in polar coordinates.

Assuming that is bounded, let us solve the following IVBVP in a circular membrane with radius 1:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

The Vibrating Drumhead

Suppose that we want to study the vibrations of a membrane of a circular drum of radius a. Then the two-dimensional wave equation [150]:

becomes:

with in polar coordinates.

Assuming that is bounded, let us solve the following IVBVP in a circular drum with radius 2:

If we truncate the infinite summation in the solution at, for example, , we get:

Now we visualize the three-dimensional plot of the solution:

We visualize the three-dimensional plot of the solution in as:

Vibrations of a Circular Membrane with Angular Periodicity

If we consider the vibrations of a circular membrane with angular periodicity, for the following initial and boundary value problem, assuming that is bounded, DSolveValue returns [151]:

For practical purposes, we truncate the infinite summation at, for example, in the preceding solution. Hence we need to compute the first 50 eigenvalues for and :

Now compute the summation using these 50 eigenvalues:

We visualize the three-dimensional plot of the solution as:

We visualize the three-dimensional plot of the solution in as:

The Potential or Laplace Equation

The Laplace equation for the steady-state temperature distribution in two dimensions is

Equation (152) also expresses the time-independent equilibrium state of a two-dimensional membrane, hence it is a substantial common part of both the heat and wave equations in two dimensions. The Laplacian of a function tells us about the deviation of the function's value at a point from the average of its values at neighboring points. For instance, a steady-state, stretched stationary rubber membrane satisfies the Laplace equation, thus the height of the membrane at any point is equal to the average height of the membrane at nearby points [153, 154].

Poisson's equation is =f, where f is a function of the space variables alone, and characterizes several physical phenomena such as [155]:

The analogous Laplace equation in three dimensions is:

Potential in a Rectangle

Dirichlet's problem in a rectangle is one of the elementary and notable problems in mathematical physics. Suppose that we would like to study the electrostatic potential in a charge-free rectangular domain [156]:

For this BVP, DSolveValue returns:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Further Potential BVP Examples in a Rectangle

In this subsection, we consider examples of BVPs with conditions on the boundary given as Dirichlet conditions, Neumann conditions or Robin conditions.

Suppose that the BCs involve both the Dirichlet and the Neumann conditions [157]:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Suppose that the BCs involve the Dirichlet, Neumann and Robin conditions [158]:

For practical purposes, we truncate the infinite summation at, for example, in the preceding solution. Hence we need to compute the first 200 eigenvalues for :

Now compute the summation using these 200 eigenvalues:

A three-dimensional plot of the truncated solution is:

Potential in Unbounded Regions

In this subsection, we consider potential problems over the upper half-plane, right half-plane, infinite strip, semi-infinite vertical strip and semi-infinite horizontal strip.

Assume that is absolutely integrable over the infinite semi-half-plane, and consider the steady-state temperature distribution on the upper half-plane. So the BVP is [159]:

for which DSolveValue returns:

A three-dimensional plot of the solution is:

Assume that is absolutely integrable over the infinite semi-half-plane, and consider the steady-state temperature distribution on the right half-plane. So the BVP is [160]:

for which DSolveValue returns:

A three-dimensional plot of the solution is:

Assume that is absolutely integrable over the infinite thin plate, and consider the steady-state temperature distribution in a thin plate whose lateral surfaces are insulated. So the BVP is [161]:

for which DSolveValue returns:

Let us compute the list of data points of the solution using NIntegrate:

Now a three-dimensional plot of the solution using ListPlot3D is:

Assume that is absolutely integrable over the semi-infinite strip, and consider the steady-state temperature distribution over this semi-infinite vertical slot. So the BVP is [162]:

for which DSolveValue returns:

A three-dimensional plot of the solution for is:

A three-dimensional plot of the solution for is:

Assume that is absolutely integrable over the semi-infinite strip, and consider the steady-state temperature distribution over this semi-infinite horizontal slot. So the BVP is [163]:

for which DSolveValue returns:

A three-dimensional plot of the solution for is:

A three-dimensional plot of the solution for is:

Inhomogeneous Boundary Conditions

Let us consider a problem with sources acting at the boundary [164]:

For this BVP with inhomogeneous BCs, DSolveValue returns:

If we truncate the infinite summation at, for example, , a three-dimensional plot of the solution is:

Laplace Equation in Polar Coordinates

In this subsection, we consider potential BVPs in polar coordinates.

Suppose that we want to solve the interior Dirichlet problem describing the electrostatic potential inside the circle when the potential is given on the boundary [165]:

For this BVP, DSolveValue returns:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Suppose that we want to solve the interior Dirichlet problem between two circles, i.e. in an annulus, when the potential is given on the boundary [166]:

For this BVP, DSolveValue returns:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Assuming that the temperature at the origin is finite, suppose we study the steady-state temperature distribution in a thin plate over a cylindrical domain with insulated lateral surface and periodic BCs [167]:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Assuming that the temperature at the origin is finite, suppose we consider the steady-state temperature distribution in a thin plate over a cylindrical domain with insulated lateral surface [168]:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Suppose that we want to find the potential in a disk in response to a forcing term acting inside the disk. Hence we have an inhomogeneous Laplace PDE-BVP with inhomogeneous BCs on a disk with radius a, for which DSolveValue returns [169]:

If , and , we get:

A three-dimensional plot of the solution is:

A three-dimensional plot of the solution in the circular region is:

Laplace Equation in Spherical Coordinates

In this subsection, we consider potential BVPs in spherical coordinates.

Suppose we want to find the steady-state temperature between two spheres held at different temperatures [170]:

For this BVP, DSolveValue returns:

For , , and , we have:

Suppose that we want to find the steady-state temperature inside a sphere when the temperature (depending on θ alone) is specified on the boundary [171]:

For this BVP, DSolveValue returns:

For and , we have:

Laplace Problem in 3D

Suppose that we consider the Laplace problem in a rectangular parallelepiped [172]:

If we truncate the infinite summation at, for example, , a three-dimensional contour plot of the solution is:

Nonlinear PDEs

As the superposition principle does not hold for nonlinear equations, and hence the method of eigenfunctions and the transform methods cannot be applied, we have a number of methods to deal with nonlinear problems.

Shock Waves and Rarefaction Fan Solutions

The IVP in which the initial data is formed of two constant states and and has a jump discontinuity at is called the Riemann problem. Now let us consider Burgers' equation with the Riemann problem [173, 174]:

where the IC is given as:

For , the characteristics intersect, hence a smooth solution cannot be constructed. Instead, we can construct a weak solution of the form of a shock wave moving at a speed determined by the Rankine–Huguenot condition. For , we can construct a weak solution of the form of a rarefaction fan (transonic rarefaction) from the simple-wave principle:

So when u(x,0)=, we have a shock wave solution:

for which a three-dimensional plot of the solution is:

Whereas when u(x,0)=, we have a rarefaction fan solution:

for which a three-dimensional plot of the solution is:

Consider the generalized viscous Burgers' equation [175] from the theory of turbulence:

where the IC is given as:

For this IVBVP, DSolveValue returns:

The solution is smooth for any positive value of ϵ:

The solution develops a shock discontinuity in the limit when ϵ approaches 0:

Charpit’s Method

Suppose that we have first-order partial differential order for , of the form [176]:

where and . Our goal is to find another equation that has the same solution as equation (177), of the form , where is a constant. Then the idea is combining these two equations and solving for p and q, which is followed by solving the Pfaffian differential equation:

for u. The requirements that and should be such that their solutions are the same and (178) is integrable lead to the compatibility condition:

The condition (179) leads to the characteristic equations (Charpit equations):

If we can find any integral of (180), involving p or q or both, we can take this integral as an additional differential equation, namely , and we can solve:

for p and q to make (181) integrable. Integrating (182) gives the solution to (183).

Note that solving the Pfaffian differential equation (184) leads to a complete integral of (185). The function CompleteIntegral returns directly the complete integrals of first order PDEs.

Suppose we have the following nonlinear differential equation:

Then from (186), we have

Using the equalities in (187), we get:

Solving (188) for u with our findings from (189) gives us the solution that we are after:

DSolve gives the same solution for this PDE, agreeing with (190):

CompleteIntegral gives the same complete integral:

Jacobi’s Method

Suppose that we have first-order partial differential order for , of the form [191]:

where , and . Our goal is to find two other equations that have the same solution as equation (192), of the form and . Then the idea is combining these three equations and solving for , and , which is followed by solving the Pfaffian differential equation:

for u. The requirement that and should be such that their solutions are the same as the solution to (193) leads to the compatibility condition:

where is the Poisson bracket.

The conditions (194) lead to the characteristic equations:

If we can find integrals of (195), namely and , and then verify that , we can solve:

for , and , which will make (196) integrable. Integrating (197) gives the solution to (198).

The method is applicable to first-order PDEs with three or more independent variables.

Similar to the Charpit's method, solving the Pfaffian differential equation (199) leads to a complete integral of (200), which can be computed directly by CompleteIntegral.

Suppose we have the following nonlinear differential equation:

Hence from (201) we have:

Note that . Solving (202) and (203) for , and gives:

Using our findings from (204), the Pfaffian differential equation (205) becomes:

Integrating (206) gives:

DSolve gives the same solution for this PDE, agreeing with (207):

CompleteIntegral gives the same complete integral:

Traveling Wave Solutions

In a traveling frame of reference, , by letting the new independent variable, say, , one can transform a polynomial PDE in into an ODE in the new independent variable . Since the derivative of is polynomial in , i.e. , all derivatives of are polynomials of . Therefore, using the chain rule, the polynomial PDE is transformed into an ODE in , which has polynomial coefficients in . One then can seek polynomial solutions of the ODE, and thus generate a subset of the set of all solutions. The argument here could be generalized to cover solutions in other hyperbolic trigonometric functions [208].

The ubiquitous KdV equation:

models shallow water waves in a canal and ion-acoustic waves in plasmas [209]:

Suppose an initial and a boundary condition are present:

Here is a visualization of the plot of the solution:

The classical Boussinesq equation,

with real parameter , was proposed by Boussinesq to describe surface water waves whose horizontal scale is much larger than the depth of the water [210]:

Suppose for an initial and a boundary condition are present:

Here is a visualization of the plot of the solution:

The Kadomtsev–Petviashvili (KP) equation is a nonlinear partial differential equation in two spatial coordinates and one temporal coordinate modeling the evolution of nonlinear long waves of small amplitude with slow dependence on the transverse coordinate [211]:

The case is known as the KPII equation and models, for instance, water waves with small surface tension. The case is known as the KPI equation and may be used to model waves in thin films with high surface tension.

The KP equation is regarded as a universal integrable equation in two spatial dimensions in the same way that the KdV equation is regarded as a universal integrable equation in one spatial dimension. The KP equation is often written as a system in compatibility form of its Lax pair:

for which DSolve returns:

Suppose for KPII (i.e. ), we have an IC and BCs:

Here is a visualization of the three-dimensional plot of the solution:

The KdV–ZK equation [212],

models ion-acoustic waves in magnetized multi-component plasmas, including negative ions:

Suppose for , an IC and BCs are present:

Applications

In this section, we consider applications of PDE-IVBVPs to problems from physics, engineering, finance, etc.

The Vibrating Beam (Fourth-Order PDE)

As opposed to the transverse vibrations of a violin string, the transverse vibrations of a thin beam show resistance to bending, and this resistance changes the wave equation to the fourth-order beam equation:

where is the ratio of rigidity constant to the linear density of the beam [213].

If we consider the small vibrations of a thin beam whose ends are held stationary, but the slopes at the endpoints can move, then the bending moments of the beam, represented by , should be zero at the endpoints. Then we should have the following initial and boundary conditions [214]:

As an example, suppose that we have a beam with length 1 and ; then the following PDE, BCs and ICs make up the IVBVP:

for which DSolveValue returns:

A three-dimensional plot of the solution is:

The Black–Scholes Equation

All European options satisfy the Black–Scholes PDE when the underlying asset price is assumed to be risk-neutral. If is the value of the option with underlying asset and elapsed time , risk-free interest rate , dividend yield and volatility or annualized standard deviation of the asset price , then the equation modeling all European options is [215]:

If is the strike and is the time of maturity, the IC for the put option is given as [216]:

For the put option problem, DSolveValue returns:

For the put option, the BCs are implicitly taken as:

Compute the value of an option for typical values of the parameters:

Compare with the value given by FinancialDerivative:

The IC for the call option is given as [217]:

For the call option problem, DSolveValue returns:

For the call option, the BCs are implicitly taken as:

Compute the value of an option for typical values of the parameters:

Compare with the value given by FinancialDerivative:

The Schrödinger Equation

In quantum mechanics, motion of a free particle with mass m is described by the time-dependent free Schrödinger equation:

where is the Planck constant divided by and u is the wavefunction describing the motion of the particle [218].

Then for the following IVBVP, DSolveValue returns:

Assuming , the three-dimensional plot of the solution follows:

Spherical Harmonics

Suppose that we want to find the potential inside a sphere with radius 3 when the potential, , is given on the boundary. For this we have to solve the Laplace PDE-BVP in a sphere with no azimuthal symmetry [219]:

Let us define the truncated solution:

Notice that truncating the summation at larger values of also gives the same solution:

Here is a visualization of the spherical three-dimensional plot of the solution:

Systems of PDEs

Many physical systems cannot be described by a single PDE but, in fact, are modeled by a system of PDEs. In these equations, unknown functions like pressure, density and temperature and their partial derivatives are described by physical laws, and we aim to find all of these functions at the same time [220].

Yet another reason for studying systems of PDEs is that it is often possible to write a single higher-order PDE as a system of first-order partial differential equations [221].

For the following IVBVP for a system of PDEs with constant coefficients, DSolve returns:

With , and , we have:

Here is the three-dimensional plot of :

Here is the three-dimensional plot of :

For the following IVBVP for an inhomogeneous system of PDEs with variable coefficients, DSolveValue returns:

Here is the three-dimensional plot of :

Here is the three-dimensional plot of :

Modeling

Heat transfer in thermal engineering is interested in modeling of transmission of energy. Let us solve a slightly modified version of the basic heat transfer example:

The model domain Ω of width and height is a ceramic strip that is embedded in a high-thermal-conductive material. The side boundaries of the strip are maintained at a constant temperature . The top surface of the strip is losing heat via thermal convection, and thermal radiation to the ambient environment at is omitted. The bottom boundary, however, is assumed to be thermally insulated.

Our goal is to find the steady-state temperature distribution of the ceramic strip.

Set up a rectangular domain with a width of and a height of .

The thermal conductivity k, heat transfer coefficient h, density ρ, heat capacity and emissivity ε of the ceramic trip are given by:

Define the model parameters:

Set up temperature surface boundary conditions at the left and right boundaries:

Set up the convective boundary condition on the top surface:

A default thermally insulated boundary condition is implicitly applied on the remaining bottom boundary:

Solve the heat transfer PDE model with NDSolveValue and DSolveValue:

Visualize the steady-state temperature distribution:

Summary

An overview of PDEs, IVBVPs and methods available in DSolve and DSolveValue to solve these problems are given in this monograph.

References

1.  Strauss, W. Partial Differential Equations: An Introduction. (2nd ed.) John Wiley & Sons, Inc., 2008.

2.  Farlow, S. J. Partial Differential Equations for Scientists and Engineers. Dover Publications, Inc., 1993.

3.  Powers, D. L. Boundary Value Problems and Partial Differential Equations. (6th ed.), 2010.

4.  Trim, D. W. Introduction to Complex Analysis and Its Application. (1st ed.) PWS Pub. Co., 1996.

5.  Evans, L. C. Partial Differential Equations. (2nd ed.) American Mathematical Society, 2010.

6.  Zwillinger, D. Handbook of Differential Equations. (3rd ed.) Academic Press, 1997.

7.  Baldwin, D., Göktaş Ü., Hereman W., Hong L., Martino R. S. and Miller J. C. "Symbolic Computation of Exact Solutions Expressible in Hyperbolic and Elliptic Functions for Nonlinear PDEs." Journal of Symbolic Computation 37, no. 6 (2004): 669–705.

8.  Kelly, M. "Evaluation of Financial Options Using Radial Basis Functions in Mathematica." The Mathematica Journal 11, no. 3 (2010): 333–257.

9.  Deconinck, B., Drogdon T. and Vasan V. "The Method of Fokas for Solving Linear Partial Differential Equations." SIAM Review 56 (2014): 159–186.

10.  Stavroulakis I. P. and Tersian S. A. Partial Differential Equations—An Introduction with Mathematica and Maple. (2nd ed.) World Scientific Pub. Co. Re. Ltd., 1999.