# Introduction to Advanced Numerical Differential Equation Solving in *Mathematica*

## Overview

The *Mathematica* function NDSolve is a general numerical differential equation solver. It can handle a wide range of *ordinary differential equations* (ODEs) as well as some *partial differential equations* (PDEs). In a system of ordinary differential equations there can be any number of unknown functions , but all of these functions must depend on a single "independent variable" , which is the same for each function. Partial differential equations involve two or more independent variables. NDSolve can also solve some *differential-algebraic equations* (DAEs), which are typically a mix of differential and algebraic equations.

NDSolve[{eqn_{1},eqn_{2},...},u,{t,t_{min},t_{max}}] | find a numerical solution for the function u with t in the range to |

NDSolve[{eqn_{1},eqn_{2},...},{u_{1},u_{2},...},{t,t_{min},t_{max}}] | find numerical solutions for several functions |

Finding numerical solutions to ordinary differential equations.

NDSolve represents solutions for the functions as InterpolatingFunction objects. The InterpolatingFunction objects provide approximations to the over the range of values to for the independent variable .

In general, NDSolve finds solutions iteratively. It starts at a particular value of , then takes a sequence of steps, trying eventually to cover the whole range to .

In order to get started, NDSolve has to be given appropriate initial or boundary conditions for the and their derivatives. These conditions specify values for , and perhaps derivatives , at particular points . When there is only one at which conditions are given, the equations and initial conditions are collectively referred to as an initial value problem. A boundary value occurs when there are multiple points . NDSolve can solve nearly all initial value problems that can symbolically be put in normal form (i.e. are solvable for the highest derivative order), but only linear boundary value problems.

In[1]:= |

Out[1]= |

When you use NDSolve, the initial or boundary conditions you give must be sufficient to determine the solutions for the completely. When you use DSolve to find symbolic solutions to differential equations, you may specify fewer conditions. The reason is that DSolve automatically inserts arbitrary symbolic constants C[i] to represent degrees of freedom associated with initial conditions that you have not specified explicitly. Since NDSolve must give a numerical solution, it cannot represent these kinds of additional degrees of freedom. As a result, you must explicitly give all the initial or boundary conditions that are needed to determine the solution.

In a typical case, if you have differential equations with up to derivatives, then you need to either give initial conditions for up to derivatives, or give boundary conditions at points.

In[2]:= |

Out[2]= |

In[3]:= |

Out[3]= |

In[4]:= |

Out[4]= |

You can use NDSolve to solve systems of coupled differential equations as long as each variable has the appropriate number of conditions.

In[5]:= |

Out[5]= |

In[6]:= |

Out[6]= |

You can give initial conditions as equations of any kind. If these equations have multiple solutions, NDSolve will generate multiple solutions.

In[7]:= |

Out[7]= |

NDSolve was not able to find the solution for y'[x]=-Sqrt[y[x]^3], because of problems with the branch cut in the square root function.

In[8]:= |

Out[8]= |

NDSolve can solve a mixed system of differential and algebraic equations, referred to as *differential-algebraic equations* (DAEs). In fact, the example given is a sort of DAE, where the equations are not expressed explicitly in terms of the derivatives. Typically, however, in DAEs, you are not able to solve for the derivatives at all and the problem must be solved using a different method entirely.

In[9]:= |

Out[9]= |

Note that while both of the equations have derivative terms, the variable appears without any derivatives, so NDSolve issues a warning message. When the usual substitution to convert to first-order equations is made, one of the equations does indeed become effectively algebraic.

Also, since only appears algebraically, it is not necessary to give an initial condition to determine its values. Finding initial conditions that are consistent with DAEs can, in fact, be quite difficult. The tutorial "Numerical Solution of Differential-Algebraic Equations" has more information.

In[10]:= |

Out[10]= |

From the plot, you can see that the derivative of is tending to vary arbitrarily fast. Even though it does not explicitly appear in the equations, this condition means that the solver cannot continue further.

Unknown functions in differential equations do not necessarily have to be represented by single symbols. If you have a large number of unknown functions, for example, you will often find it more convenient to give the functions names like or .

In[11]:= |

Out[16]= |

This actually computes an approximate solution of the heat equation for a rod with constant temperatures at either end of the rod. (For more accurate solutions, you can increase .)

In[17]:= |

Out[17]= |

An unknown function can also be specified to have a vector (or matrix) value. The dimensionality of an unknown function is taken from its initial condition. You can mix scalar and vector unknown functions as long as the equations have consistent dimensionality according to the rules of *Mathematica* arithmetic. The InterpolatingFunction result will give values with the same dimensionality as the unknown function. Using nonscalar variables is very convenient when a system of differential equations is governed by a process that may be difficult or inefficient to express symbolically.

In[18]:= |

Out[19]= |

NDSolve is able to solve some partial differential equations directly when you specify more independent variables.

NDSolve[{eqn_{1},eqn_{2},...},u,{t,t_{min},t_{max}},{x,x_{min},x_{max}},...] | |

solve a system of partial differential equations for a function with t in the range to and x in the range to , ... | |

NDSolve[{eqn_{1},eqn_{2},...},{u_{1},u_{2},...},{t,t_{min},t_{max}},{x,x_{min},x_{max}},...] | |

solve a system of partial differential equations for several functions |

Finding numerical solutions to partial differential equations.

In[20]:= |

Out[20]= |

In[21]:= |

Out[21]= |

NDSolve currently uses the numerical method of lines to compute solutions to partial differential equations. The method is restricted to problems that can be posed with an initial condition in at least one independent variable. For example, the method cannot solve elliptic PDEs such as Laplace's equation because these require boundary values. For the problems it does solve, the method of lines is quite general, handling systems of PDEs or nonlinearity well, and often quite fast. Details of the method are given in "Numerical Solution of Partial Differential Equations".

In[22]:= |

Out[22]= |

In[23]:= |

Out[23]= |

As mentioned earlier, NDSolve works by taking a sequence of steps in the independent variable t. NDSolve uses an adaptive procedure to determine the size of these steps. In general, if the solution appears to be varying rapidly in a particular region, then NDSolve will reduce the step size to be able to better track the solution.

NDSolve allows you to specify the precision or accuracy of result you want. In general, NDSolve makes the steps it takes smaller and smaller until the solution reached satisfies either the AccuracyGoal *or* the PrecisionGoal you give. The setting for AccuracyGoal effectively determines the absolute error to allow in the solution, while the setting for PrecisionGoal determines the relative error. If you need to track a solution whose value comes close to zero, then you will typically need to increase the setting for AccuracyGoal. By setting AccuracyGoal->Infinity, you tell NDSolve to use PrecisionGoal only. Generally, AccuracyGoal and PrecisionGoal are used to control the error local to a particular time step. For some differential equations, this error can accumulate, so it is possible that the precision or accuracy of the result at the end of the time interval may be much less than what you might expect from the settings of AccuracyGoal and PrecisionGoal.

NDSolve uses the setting you give for WorkingPrecision to determine the precision to use in its internal computations. If you specify large values for AccuracyGoal or PrecisionGoal, then you typically need to give a somewhat larger value for WorkingPrecision. With the default setting of Automatic, both AccuracyGoal and PrecisionGoal are equal to half of the setting for WorkingPrecision.

NDSolve uses error estimates for determining whether it is meeting the specified tolerances. When working with systems of equations, it uses the setting of the option NormFunction->f to combine errors in different components. The norm is scaled in terms of the tolerances, given so that NDSolve tries to take steps such that

where is the component of the error and is the component of the current solution.

In[24]:= |

Out[24]= |

In[25]:= |

Out[25]= |

Through its adaptive procedure, NDSolve is able to solve "stiff" differential equations in which there are several components varying with at extremely different rates.

NDSolve follows the general procedure of reducing step size until it tracks solutions accurately. There is a problem, however, when the true solution has a singularity. In this case, NDSolve might go on reducing the step size forever, and never terminate. To avoid this problem, the option MaxSteps specifies the maximum number of steps that NDSolve will ever take in attempting to find a solution. For ordinary differential equations, the default setting is MaxSteps->10000.

In[26]:= |

Out[26]= |

In[27]:= |

Out[27]= |

The default setting MaxSteps->10000 should be sufficient for most equations with smooth solutions. When solutions have a complicated structure, however, you may sometimes have to choose larger settings for MaxSteps. With the setting MaxSteps->Infinity there is no upper limit on the number of steps used.

NDSolve has several different methods built in for computing solutions as well as a mechanism for adding additional methods. With the default setting Method->Automatic, NDSolve will choose a method which should be appropriate for the differential equations. For example, if the equations have stiffness, implicit methods will be used as needed, or if the equations make a DAE, a special DAE method will be used. In general, it is not possible to determine the nature of solutions to differential equations without actually solving them: thus, the default Automatic methods are good for solving as wide variety of problems, but the one chosen may not be the best one available for your particular problem. Also, you may want to choose methods, such as symplectic integrators, which preserve certain properties of the solution.

Choosing an appropriate method for a particular system can be quite difficult. To complicate it further, many methods have their own settings, which can greatly affect solution efficiency and accuracy. Much of this documentation consists of descriptions of methods to give you an idea of when they should be used and how to adjust them to solve particular problems. Furthermore, NDSolve has a mechanism that allows you to define your own methods and still have the equations and results processed by NDSolve just as for the built-in methods.

When NDSolve computes a solution, there are typically three phases. First, the equations are processed, usually into a function that represents the right-hand side of the equations in normal form. Next, the function is used to iterate the solution from the initial conditions. Finally, data saved during the iteration procedure is processed into one or more InterpolatingFunction objects. Using functions in the context, you can run these steps separately and, more importantly, have more control over the iteration process. The steps are tied by an object, which keeps all of the data necessary for solving the differential equations.