Wolfram Computation Meets Knowledge

New in the Wolfram Language: Symbolic PDEs

PDEs word cloud

Partial differential equations (PDEs) play a vital role in mathematics and its applications. They can be used to model real-world phenomena such as the vibrations of a stretched string, the flow of heat in a bar, or the change in values of financial options. My aim in writing this post is to give you a brief glimpse into the fascinating world of PDEs using the improvements for boundary value problems in DSolve and the new DEigensystem function in Version 10.3 of the Wolfram Language.

The history of PDEs goes back to the works of famous eighteenth-century mathematicians such as Euler, d’Alembert, and Laplace, but the development of this field has continued unabated during the last three centuries. I have, therefore, chosen examples of both classical as well as modern PDEs in order to give you a taste of this vast and beautiful subject.

Let us begin by considering the vibrations of a stretched string of length π that is fixed at both ends. The vibrations of the string can be modeled by the one-dimensional wave equation, which is given below. Here, u(x,t) is the vertical displacement of a point at the location x on the string, at time t:

Vibrations of the string can be modeled by the one-dimensional wave equation

Next, we specify a boundary condition to indicate that the ends of the string remain fixed during the vibrations:

Specifying a boundary condition

Finally, we prescribe the displacement and speed at different points on the string at time t=0 as initial conditions for the motion of the string:

Prescribing initial conditions for the motion of the string

We can now use DSolve to solve the initial boundary value problem for the wave equation:

Using DSolve to solve the initial boundary value problem for the wave equation

As seen above, the solution is an infinite sum of trigonometric functions. The sum is returned in an Inactive form since each individual term in the expansion has a physical interpretation and, typically, a small number of terms will provide a useful approximation to the exact solution. For example, we can extract four terms from the sum to obtain an approximate solution asol(x,t):

Extracting four terms from the sum

Each term in the sum represents a standing wave, which can be visualized as follows:

Visualizing a standing wave

These standing waves combine in a magical way to produce the smooth vibrations of the string, as shown in the following animation:

Animation of smooth vibrations of the string

The wave equation belongs to the class of linear hyperbolic PDEs, which govern the propagation of signals with finite speeds. This PDE arose as a convenient way to model the vibrations of strings and other deformable bodies, but it plays an even more important role in modern physics and engineering since it governs the propagation of light and other electromagnetic waves.

Let us now model the flow of heat in a bar of length 1 that is insulated at both ends using the heat equation, which is given as follows:

Modeling the flow of heat using the heat equation

Since the bar is insulated at both ends, no heat flows through the ends, which translates into a boundary condition at the two ends x=0 and x=1:

Boundary condition at the two ends x=0 and x=1

Next, we must prescribe the initial temperature at different points in the bar. In this example, we choose the linear function given below. The initial temperature is 20 at the left end (x=0) and 100 at the right end (x=1):

initialcondition = u[x, 0] = 20 + 80 x

Finally, we solve the heat equation subject to these conditions:

Solving the heat equation

As in the case of the wave equation, we can extract a few terms from the Inactive sum to obtain an approximate solution, as shown below:

Extracting a few terms for the Inactive sum to obtain an approximate solution

The first term in the approximate solution, 60, is the average of the initial temperatures at the left and the right ends of the bar, and is also the steady-state temperature for this bar. As the plot of temperature along the length of the bar below shows, the temperature of the bar quickly evolves to the steady-state value of 60 degrees:

Plot of temperature along the length of the bar

The heat equation belongs to the class of linear parabolic PDEs, which govern diffusion processes. This rather simple-looking equation is ubiquitous and makes a surprising appearance in a wide variety of applications. We will see two examples of this phenomenon later in the post.

We now turn to the Laplace equation, which is used to model the steady state of systems, i.e. the behavior after any time-dependent transients have died out. In two dimensions, this equation is given as follows:

Laplace equation in two dimensions

Assume that the coordinates x and y are restricted to the rectangular region Ω given below:

Omega = Rectangle [{0,0}, {1,2}]

The classical Dirichlet problem is to find a function u(x,y) that satisfies the Laplace equation inside Ω along with a given DirichletCondition that specifies values on the boundary of Ω, as illustrated below:

Classic Dirichlet problem

The Dirichlet problem can be solved using the elegant region notation for DSolve:

The Dirichlet problem can be solved using the elegant region notation for DSolve

As in earlier examples, we can extract a fixed number, say 100, of terms from the Inactive sum and visualize the solution, as shown below:

Visualizing the solutions from the Inactive sum

Notice that the solution u(x,y) of the Dirichlet problem is smooth within Ω despite the fact that the boundary conditions have sharp features. Also, u(x,y) attains its maximum and minimum values on the boundary, while there is a saddle point at the center of the rectangle. These features are typical of linear elliptic equations, the class of PDEs to which the Laplace equation belongs.

The wave, heat, and Laplace equations are the most famous examples of classical PDEs. We will now consider three representative examples of modern PDEs, beginning with Burgers’ equation for viscous fluid flow, which can be stated as follows:

Burgers' equation for viscous fluid flow

This nonlinear PDE was introduced by J. M. Burgers around 1940 as a simple model for turbulent flows (the parameter ϵ in the equation represents the viscosity of the fluid). However, a decade later, E. Hopf and J. D. Cole showed that Burgers’ equation can be transformed into the heat equation, thereby indicating that this equation cannot exhibit chaotic behavior. The Hopf–Cole transformation allows us to solve Burgers’ equation in closed form for an initial condition, such as the one given below:

Using Hopf--Cole transformation to solve Burgers' equation in closed form for an initial condition

For this example, we use DSolveValue, which returns only an expression for the solution. The error function (Erf) terms in the formula below arise from the solution of the corresponding initial value problem for the heat equation:

Using DSolveValue

The following plot is showing the evolution over time of a hypothetical fluid velocity field in one dimension. The solution is smooth for a positive value of ϵ, although the initial condition is a piecewise function:

Plot is showing the evolution over time of a hypothetical fluid velocity field in one dimension

As seen below, the solution develops a shock discontinuity in the limit when the viscosity ϵ approaches 0. These shock solutions are a well-known feature of the inviscid (ϵ = 0) Burgers’ equation:

Shock discontinuity in the limit when the viscosity epsilon approaches 0

As a second example of a modern PDE, let us consider the Black–Scholes PDE from finance. This equation was introduced by Fischer Black and Myron Scholes in 1973 as a model for the prices of European-style financial options, and can be stated as follows:

Black--Scholes PDE from finance

Here, c is the price of the option as a function of stock price s and time t, r is the risk-free interest rate, and σ is the volatility of the stock.

In their landmark paper (which has been cited more than 28,000 times), Black and Scholes noted that their equation can be reduced to the heat equation by a transformation of variables. This dramatic simplification leads to the celebrated Black–Scholes formula for European call options, under a terminal condition based on the “strike price” k of the asset at time t=T:

Black--Scholes formula for European call options

Armed with this formula, we can compute the value of a financial option for typical values of the parameters:

Computing the value of a financial option for typical values of the parameters

The answer agrees with the value given by the built-in FinancialDerivative function:

Verifying the answer with the FinancialDerivative function

As our third example of a modern PDE, we will solve the free Schrödinger equation for an electron that is constrained to move in a one-dimensional box of length d, with a suitable initial condition. The equation and conditions may be stated as follows:

Schrödinger equation and conditions

This example has an elementary solution that takes imaginary values due to the presence of I in the Schrödinger equation:

Example that has an elementary solution that takes imaginary values due to the presence of I

The probability density function for the electron, ρ=Ψ Ψ, can be computed using suitable values for the parameters in the problem as follows:

Computing the probability density function for the electron

We can create a movie showing the evolution of the probability density over time, which reveals that the “center” of the electron moves from side to side in the box:

Evolution of the probability density over time

Eigenvalues and eigenfunctions play an important role in the solution of Schrödinger’s equation and other PDEs. In particular, they provide the building blocks for the infinite sum solutions of the wave and heat equations that we saw earlier in the post. Hence, as our final example, let us consider the problem of finding the nine smallest eigenvalues and eigenfunctions for the Laplacian operator with a homogeneous (zero) Dirichlet condition over a ball in three dimensions. That is, we are interested in the nine smallest values of λ and the corresponding functions ϕ that satisfy 𝓛ϕ = λ ϕ, with the following definitions:

Finding the nine smallest eigenvalues and eigenfunctions

The new DEigensystem function in Version 10.3 allows us to compute the required eigenvalues and eigenfunctions as follows:

Version 10.3 allowing to compute the required eigenvalues and eigenfunctions

The eigenvalues for this problem are expressed in terms of BesselJZero. For example:

Eigenvalues are expressed in terms of BesselJZero

The eigenfunctions can be visualized using DensityPlot3D, which returns the beautiful plots that are shown below:

Eigenfunctions visualized using DensityPlot3D

Partial differential equations are an essential tool in many branches of science, engineering, statistics, and finance. At a more fundamental level, they provide precise mathematical formulations for some of the deepest and most subtle questions about our universe, such as whether naked singularities can really exist. In my experience, the study of PDEs rewards one with a rare combination of practical insights and intellectual satisfaction.

I invite you to look at the documentation for DSolve, NDSolve, DEigensystem, NDEigensystem, and the finite element method to learn more about the various approaches for solving PDEs in the Wolfram Language.

Symbolic PDEs is supported in Version 10.3 of the Wolfram Language and Mathematica, and is rolling out soon in all other Wolfram products.

Download this post as a Computable Document Format (CDF) file.

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.

8 comments

  1. Those are great examples and visualizations of PDEs. Eigenvalue analyses show up in all sorts of stretched strings — including the failure of the Tacoma Narrows Bridge: see the great example https://www.ma.utexas.edu/mp_arc/c/13/13-25.pdf .
    I don’t understand the historic commentary about chaotic behavior and PDEs. A quick scan of the literature shows that PDEs are used to model chaotic behavior in all sorts of applications.

    Reply
    • The author was referring specifically to Burger’s equation, which is a *simple* model of turbulence. Not all nonlinear PDEs exhibit chaos.

      Reply
    • Thank you for the Tacoma Narrows Bridge reference.

      As noted by you, nonlinear PDEs are certainly used to model chaotic behavior. However, since the initial value problem for Burgers’ equation can be solved in closed form using the Hopf-Cole transformation, this nonlinear PDE does not provide a good model for studying chaotic behavior. I hope that this clarifies the comment made in the post.

      Reply
  2. In addition to learning about the new PDE capabilities I also learned, quite to my surprise, that I could pattern match the Infinity symbol, i.e. Infinity->6 for example..

    Reply
  3. Hello Devendra,

    How we can solve heat equation not in disk. As example in ring or outer zone of disk. These examples have difficulties in solving because of definition of eigenvalues is very hard.

    Reply
    • Hello Andrey,

      Thank you for the comment. We expect to add support for solving the heat equation in a thin circular ring and other regions using DSolve in a future release, along with related functionality such as solving Sturm-Liouville problems with periodic boundary conditions.

      Reply