This article describes how to model diffusion using NDSolve, and then compares that to constructing your own methods using procedural, functional, rule-based, and modular programming. While based on the diffusion equation, these techniques can be applied to any partial differential equation.

Introduction

Using the built-in Mathematica command NDSolve to solve partial differential equations is very simple to do, but it can hide what is really going on. You can always consult references about using Mathematica for differential equations [1, 2]. Exploring various programming methods also lets you create your own procedures that can be incorporated into NDSolve [1]. I chose the diffusion equation as the main example because there is so much material available for it and because of its high level of interest [3, 4, 5]. In this article I am using Mathematica 8.

Begin with a model of diffusion, in this case, the diffusion equation. The continuity equation is

where is the density, is time, and is the component of velocity. We can assume the rate of diffusion is constant and write the diffusion equation,

(1)

For complete generality, we can write

where is a positive-definite matrix of diffusion coefficients that may or may not depend on position and time; this is the diffusion tensor.

Building Up Your Program Idea in Small Steps

First, try built-in functions to see if they can solve your problem right away. The plan here is to use DSolve and then, when it likely fails to solve the PDE, try NDSolve. We are solving for a ring, so we will use only one spatial dimension. We begin by writing our diffusion equation.

Now we try DSolve.

DSolve is unable to reach a solution to the naively formulated problem. We next provide a set of initial conditions.

We incorporate the initial value and try again.

Again it fails; we could play with the initial conditions until we found a set that worked, or we could try NDSolve. Let us say that we want to investigate the diffusion in a tube of gas. We consider this tube to be one dimensional to keep things simple. We assume our ring has a fixed length 1. We then “unwrap” the ring, making it a line, which lets us plot results. We will solve the equation for 10 units of time. We also specify the rate of diffusion.

Let us try changing the initial density value to 0.

We put all of these elements together into NDSolve and get a solution in the form of an interpolating function.

We can plot this function.

Figure 1. Solution to the diffusion equation with initial density of 0 in empty space.

Let us try another initial value, say a sinusoidal density wave.

Now we try a solution.

We can plot this.

Figure 2. Solution to the diffusion equation with initial density based on a sine function. We have not determined the rate of diffusion.

This does not seem realistic, as the density drops to zero immediately. Now we need to determine what is, say 1.

Here is the plot.

Figure 3. Solution to the diffusion equation with sinusoidal boundary conditions.

This seems more realistic than Figure 2, but the boundary conditions do not match up. If we think of this as a circle (wrapping the line to form a ring), we suddenly get a discontinuity when we go from to . We can repair our ring solution by using periodic boundary conditions.

We are forcing an initial uniform density of 1, despite the sinusoidal boundary conditions in time, which would give us values of 0 density at the origin, and we do not want that. We get the solution, along with a warning message about inconsistent boundary and initial value conditions. Please ignore that in this case.

We can plot this.

Figure 4. A fairly good model of diffusion on a ring with sinusoidal boundary conditions.

This looks reasonable. The boundary conditions are now periodic, as on a ring. How do we write a program to do this?

Procedural Programming of Diffusion

Those of you who have programmed in languages like FORTRAN, C++, or Java are familiar with the idea of procedural programming. Procedures are executed in looping structures. We implement a finite-difference scheme to solve our equation. NDSolve is able to use finite differences as a specific method [1].

Begin by converting the diffusion equation to its finite-difference equivalent,

(2)

Here a ring is divided into cells for time steps. Specify a given cell’s density by considering the cell for the time step. So extends from 0 to , and time extends from 0 to . For any specific time step,

The stability condition [2] is

(3)

If this ratio is greater than 1/2, the solution becomes unstable in time. Write

(4)

If 0 is some set of initial conditions for , in the next time step we can use this equation to get a set of values for . We can then relabel these as the new 0 values and continue the calculation. We have to establish the initial conditions and use the procedural iterator to produce this result.

Choose this for positions 0 to 6.

We then develop the current time step as the initial conditions. This specifies the part of the list considered, using .

Here are the periodic boundary conditions.

This completes the initialization step, establishing the initial conditions.

We now construct a single time step by rewriting to include the boundary conditions, using .

We test this.

Define a time step, again using Table to generate the list.

We check this.

We now construct the solution in time, using Table.

This completes the model development.

Now we need to display the results. I make a test called run1.

Here is the plot of these results.

Figure 5. A jumbled mess of a solution for diffusion, not much like Figure 4.

This distorts and gets noisy very fast; of course we violated the stability condition with our choice of , so let us try it again. We know from (3) that . We need to reinitialize 0.

Then we construct run using the correct value for .

This produces the result.

Figure 6. A better solution that looks like Figure 4.

That looks better. This looks like the result we got from NDSolve, with some extra noise. We could use a finer grid and a better difference scheme if we wanted to get more precise and accurate. That is, however, the subject for another article.

Functional Programming

Another style of programming uses as the application of to . So we start by establishing a function in the traditional way, using . Here we use the notation # to represent a formal argument.

We could write this in short form without the word Function, but then we need an ampersand following the function statement so that Mathematica knows that the function statement is complete.

We write #1 and #2 for two arguments.

Or again we can write this in short form.

We can apply a function to each element of a list by the use of the command.

We could also use the short form .

This is equivalent.

Say we want to program the diffusion equation using functional programming. We begin by defining the list of cells in the ring. We use to produce a list from 1 to .

We decide we want to go from 1 to 6.

We now define the initial conditions.

For our choice of initial conditions, we have 1 everywhere, so we just say to divide by itself at each location.

Now we specify the diffusion equation. Since we have periodic boundary conditions and a list, we can just rotate the list right and left.

We can test this. First, clear the procedural definition of from the previous section.

We now need to impose our boundary conditions. Here we are replacing the end parts with the boundary conditions, using .

We will test these.

So we construct a single step.

Let us try it.

We can construct a table to produce the result.

We can try this.

We can plot this result.

Figure 7. Another solution that looks like Figure 4.

We have obtained the same result as with procedural programming.

Rule-Based Programming and Pattern-Matching

Another style of programming uses the symbolic nature of Mathematica to transform expressions. Rule-based programming uses the notion of transformation rules, which rewrite expressions based on their form. Take the expression and apply a transformation rule to it.

Here is a more complicated example.

Coming back to the diffusion equation, we begin with the initial conditions.

Then we execute it.

Then we have the diffusion equation. Here the double underscore __ following the 0 signifies a list.

For our example, define t1.

Define the boundary conditions.

Test this for the first time step.

We can test this again.

We need to include an intermediate function that specifies the values for our model for the FoldList command.

We can then produce our model runt1.

We can plot this.

Figure 8. Yet another solution that looks like Figure 4.

Scoping Constructs

A scoping construct makes symbols local so that they do not create or redefine global values. We can make local object names using the command . We can temporarily assign values to variables with the command . I will demonstrate the use of the Module command.

We can build a module for what we have been doing.

We can run this and plot the result.

Figure 9. Yet again another solution that looks like Figure 4.

Conclusion

We have seen that we can apply the canonical programming paradigms of Mathematica to the problem of diffusion. I have only scratched the surface of the diffusion problem. The methods I have presented could easily form the template for new methods for NDSolve, whose implementation is documented [1]. This flexibility in being able to introduce new methods into the existing structure of a Mathematica function is extremely powerful. Of course, the most powerful programming methods merge the paradigms.

References

[1]	M. Sofroniou and R. Knapp, Advanced Numerical Differential Equation Solving in Mathematica, Champaign, IL: Wolfram Research, Inc., 2008. www.wolfram.com/learningcenter/tutorialcollection/AdvancedNumericalDifferentialEquationSolvingInMathematica.
[2]	V. G. Ganzha and E. V. Vorozhtsov, Numerical Solutions for Partial Differential Equations—Problem Solving Using Mathematica, Boca Raton, FL: CRC Press, Inc., 1996.
[3]	R. Ghez, A Primer of Diffusion Problems, New York: John Wiley and Sons, Inc., 1988. onlinelibrary.wiley.com/book/10.1002/3527602836.
[4]	E. L. Cussler, Diffusion: Mass Transfer in Fluid Systems, 3rd ed., New York: Cambridge University Press, 2009.
[5]	J. Crank, The Mathematics of Diffusion, 2nd ed., Oxford, England: Clarendon Press, 1975.
G. E. Hrabovsky, “Diffusion Modeling,” The Mathematica Journal, 2012. dx.doi.org/doi:10.3888/tmj.14-6.

About the Author

George Hrabovsky is the president and founder of Madison Area Science and Technology, a nonprofit scientific research and education organization. He has been a Mathematica user and programmer for more than 20 years. He does research into theoretical physics, atmospheric science, and computational physics.

George E. Hrabovsky
105 Alhambra Place #2
Madison, WI 53713
george@madscitech.org