Further Examples: NDSolve
Ordinary differential equations: Basic usage
This command finds a numerical approximation to a function that is equal to its first derivative at each point x between and , and that has the value when x is .
NDSolve returns a rule to replace y by an InterpolatingFunction object.
In[1]:=
Out[1]=
Mathematica expresses the closedform solution of the following interesting differential equation as a pure function.
In[2]:=
Out[2]=
We rewrite the pure function as an ordinary algebraic expression. The solution has a singularity at .
In[3]:=
Out[3]=
Now we solve the same equation approximately, using NDSolve instead of DSolve. As we can see from the closedform solution, the warning is appropriate.
In[4]:=
Out[4]=
We can use an InterpolatingFunction object like any other function that evaluates to a number. Here are three equivalent ways to check the boundary condition.
In[5]:=
Out[5]=
In[6]:=
Out[6]=
In[7]:=
Out[7]=
It is also easy to make plots of solutions. Here is how to see a graph of an approximation to the function which is the reciprocal of twice its derivative.
Not too surprisingly, this plot looks very much like the square root function. What do you think would happen if you tried to solve in the other direction?
In[8]:=
Mathematica can handle higherorder equations. Here is a plot of the solution of a thirdorder equation.
Note that you have to specify enough initial conditions to uniquely determine the solution.
In[9]:=
Solving systems of equations works similarly. For systems of two equations a phase plot is often a good way to visualize the solution. Here is a phase plot that describes the motion of a weakly damped pendulum.
In[10]:=
In[11]:=
Ordinary differential equations: Options
This way of solving the differential equation for the cosine function doesn't work over a long interval.
In[12]:=
Out[12]=
You need to increase the setting of the MaxSteps option. As expected, at the end of the interval the solution is close to .
In[13]:=
Out[13]=
The Method option allows you to specify which method NDSolve uses to approximate the solution. The default method, Automatic, automatically switches between BDF (backward differentiation formulas) and Adams multistep methods, depending on stiffness. There are many possible settings for the Method option; different methods will be appropriate for different types of equations. A method which is generally applicable and quite effective for nonstiff equations is ExplicitRungeKutta, which uses explicit RungeKutta methods of various orders. For some problems, the RungeKutta method can find the solution using fewer steps.
For the Rössler equations, the RungeKutta method needs about half as many steps as the default method.
In[14]:=
One way to get a very precise solution of an ODE is to give a sufficiently high value for the WorkingPrecision option. Note that AccuracyGoal and PrecisionGoal default to half of the value of WorkingPrecision.
You need to be careful not to use too high a value for WorkingPrecision, because as working precision increases, not only does the time taken for each arithmetic operation and function evaluation increase, but the number of steps also typically increases exponentially. You also need to make sure that the coefficients in your problem are exact or are of at least as high a precision as the setting for WorkingPrecision.
These commands compare the solution to the Lorenz equations using different values of WorkingPrecision. They use the extrapolation stiffness switching method which is particularly effective for computing high precision solutions.
In[15]:=
In[16]:=
Out[16]=
In[17]:=
Out[17]=
In[18]:=
Out[18]=
Even though the last two solutions are given in all of the digits of the specified WorkingPrecision, they are by no means that precise. In fact, even though AccuracyGoal and PrecisionGoal default to half of the WorkingPrecision, the overall solution is not this precise because these particular equations have the property that they magnify errors (which is why they are famous!). Here is a graph using the extrapolation method displaying actual error as a function of computation time on a logarithmic scale.
The time increases exponentially with the precision.
In[19]:=
In[20]:=
Ordinary differential equations: Boundary value problems
Simple linear boundary value problems can be solved by constructing the boundary value equations appropriately. If the order of the equation is , you need to know the values of some combination of the function or its derivatives at points.
This normalized equation describes the effect of a wave incident on the right edge of an optical medium (where x is ).
In[21]:=
Here is a thirdorder equation where the function values are known at three points.
In[22]:=
In[23]:=
Differential algebraic equations
Differential algebraic equations (DAE's) are systems which are a mix of differential equations and algebraic equations to be solved simultaneously. Mathematica can solve many types of systems of differential algebraic equations.
By default, Mathematica tries to symbolically solve for the derivatives to obtain a firstorder system in the form . In the example below, two solutions are given because of the two branches of the square root function.
In[24]:=
Out[24]=
In either branch, however, the solution becomes problematic because the solution to the equation in terms of the derivative, does not have sufficient continuity near . An alternative to having NDSolve try to use symbolic solutions for the derivative is to use the option , in which case NDSolve only solves numerically to keep the equation satisfied, effectively as an algebraic condition.
In[25]:=
Out[25]=
In this case, NDSolve has solved over the entire interval, but there is only one solution because with the numerical solution, only one solution branch is followed. To solve this and other differential algebraic equations, Mathematica uses the method IDA, which was developed at the Lawrence Livermore National Laboratory.
Note that, of course, it is possible to find a solution to the equation by differentiating with respect to , which gives the very familiar equation for simple harmonic motion.
In[26]:=
Out[26]=
In general, the index of a differential algebraic system is the number of times you need to differentiate it to get a system of differential equations. The IDA method generally handles systems of index . Note, however, that by differentiating a system, you may in some cases lose some solution features and introduce numerical errors.
It is interesting to compare how well the algebraic condition is satisfied for the two successful solutions.
In[27]:=
The red is from the solution directly from the algebraic condition. It is apparent that the method is working to keep it within tolerances. On the other hand, with the differentiated solution, the algebraic equation is only satisfied indirectly, and the deviation shows drift over time.
If some of the equations you give to NDSolve are completely algebraic (they have no terms with derivatives in them), then NDSolve will solve the system using the IDA method as it does in this example from chemical kinetics.
In[28]:=
Out[28]=
In[29]:=
In[30]:=
Partial differential equations: Basic usage
For the first two examples, consider the onedimensional heat equation. With fixed boundary conditions this is a model for diffusion of heat in an insulated rod with the temperatures at the endpoints held fixed.
This command solves the heat equation with the left end () held at fixed temperature , the right end () held at fixed temperature , and an initial heat profile given by the quadratic in x.
In[31]:=
Out[31]=
Just as for ODE's, the solution is returned as a rule to replace the dependent variable with an InterpolatingFunction object. The only noticeable difference is that the InterpolatingFunction object is twodimensional. That is to say, it takes two arguments to evaluate, with the arguments in the same order as the variables in the NDSolve command.
The boundary conditions can be a linear combination of Dirichlet and Neumann type conditions and can be time dependent. For example, entering these command produces a plot of a solution with the temperature at the left edge varying sinusoidally.
In[32]:=
It is also possible to compute the solution of some systems of PDE's. These commands solve a system of mixed parabolichyperbolic type and produce separate contour plots for each of the dependent variables.
In[33]:=
There is a message indicating the error in the final solution may be large. This may occur because the number of spatial approximation points is based on the initial condition. The scaled error estimate is computed so that if it is , the AccuracyGoal and PrecisionGoal options are typically satisfied. The defaults for these are for PDE's, so even the extra factor of in error means that the error is well under 1%.
In[34]:=
In[35]:=
Periodic boundary conditions are frequently convenient for numerical solutions. You can tell NDSolve to solve with periodic boundary conditions by specifying that the values of the solution are to be equal at the left and right edges of the domain in one independent variable for all values of the other independent variable.
For example, an interesting equation with periodic solutions is the nonlinear Schrödinger equation. These commands set up a periodic initial condition (it happens to be a soliton), and computes the solution, and then produces plots showing the modulus and real and imaginary parts of the solution.
Evaluate the next cell. Then close the cell group and doubleclick on the graphic to start the animation, or select the cell bracket containing the graphics and choose Animate Selected Graphics in the Cell menu.
In[36]:=
In[37]:=
Partial differential equations: Options
Solutions for PDE's are computed using the numerical method of lines. First one variable is chosen to be a temporal variable (this must have initial conditions specified) and the others are chosen to be spatial. Next a discretization is found for the spatial variables. The result of discretizing the spatial variables is a (typically large) system of ordinary differential equations which can be solved with the many methods that NDSolve has for solving systems of ODE's. Typically the options you give directly to NDSolve will be applied solving the system of ODE's in the temporal variable. However, the following options will affect both the spatial discretization and the solution of the temporal system of ODE's: AccuracyGoal, PrecisionGoal, MaxStepFraction, MaxSteps, MaxStepSize, and StartingStepSize.
There are also options which are specific to the spatial discretization. Some of these coincide with the options listed above, allowing you to specify separate values for the spatial and temporal phases of the solution. Others are appropriate only to the spatial discretization.
This shows the possible tensor product grid options with the default settings.
In[38]:=
Out[38]=
Detailed information on how to use these options is given in the advanced documentation for NDSolve.
In the discretization stage, the default used is fourthorder finite differences. In some cases a different choice might be more efficient or give a more accurate solution. An example of this is in solving a twodimensional generalization of the sineGordon equation with periodic boundary conditions in the x and y directions.
In[39]:=
Out[39]=
This shows a plot of the solution at .
In[40]:=
For smooth solutions with periodic boundary conditions like this, it is often very efficient to use socalled pseudospectral derivatives, which are in effect a limiting case of maximal spatial difference order. Note that the timing has gone down by nearly a factor of !
In[41]:=
Out[41]=
In[42]:=
Partial differential equations: Limitations
NDSolve uses a method for solving PDE's that is called the numerical method of lines. It discretizes in one variable to make a system of ODE's. This system is then solved using the ODE methods built into NDSolve.
For the method to work, an initial function must be specified for one variable and boundary values may be specified for the other variable. The initial function is used to find the initial conditions for the system of the ODE's. Boundary and initial values may be specified on at most three sides of a rectangle.
The method has the advantage that it can solve a reasonably large class of equations. However, there are types of equations which it cannot solve.
Because elliptic problems are illposed unless boundary values are specified on all sides of a region, this method cannot find solutions of elliptic problems. A classic example is Laplace's equation. Entering this command indicates what will happen when you try to do something of this sort. Typically illposedness will appear as numerical instability, and looking at the scale on the plot produced indicates that you should heed the message seriously!
In[43]:=
Another class of problems that the method cannot currently handle are those that form singularities in the solution. The discretization is done with finite difference methods, so fronts may be incorrectly resolved or completely lost. For example, Burgers' equation is a model for some of the features of gas dynamics, including shock formation. Entering these commands produces a plot that shows the oscillations typical of the interaction between finite difference methods and fronts.
In[44]:=
In[45]:=
