Differential Equations

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

In this notebook we introduce the techniques available for solving differential equations.  The function DSolve seeks symbolic solutions while NDSolve produces numerical solutions.  The Mathematica book provides good descriptions and insightful examples of both functions and is a valuable resource but not quite a textbook.  Here we provide an introduction from the point of view of someone who has had to learn to use these functions in a more practical setting.  Although NDSolve is, in my experience, much more useful than DSolve, we begin with DSolve because it is simpler to use and many of its features carry over to NDSolve.  Then we discuss the unique features of NDSolve and some practical problems encountered in numerical work.

More information about DSolve can be found in the Mathematica book, Section 3.5.10, or in the reference guide,  DSolve.  More information about NDSolve can be found in Section 3.9.7 or NDSolve.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Needs["Utilities`Notation`"];
Needs["Graphics`"]

Symbolic solution

Basic syntax

The basic command for symbolic solution of a differential equation is DSolve.  The syntax for an ordinary differential equation with one variable is DSolve[eqn,fn,var] where eqn is a well-posed differential equation involving derivatives of the function fn with respect to the independent variable var.  The derivatives must be based upon D, not Dt.  The solution is returned in the form of a list of replacement rules for fn, with the detailed form depending on the specificity of fn.  For example, consider the equation of motion for a simple-harmonic oscillator:

sol1  =  
      DSolve[(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x,  t]

Because there exist two independent solutions to this second-order differential equation, the solution is returned in the form of a replacement rule that includes two undetermined constants, C[1] and C[2].  A more explicit solution is obtained if we specify a sufficient number of initial or boundary conditions, in the form of additional equations, that determine these constants.  In the example below, we specify the initial displacement and velocity

DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x[0] == x0,  (x^′)[0] == v0},  x,  
      t]

and obtain a complete solution free of undetermined constants.

Notice that although we must employ a dummy argument in the statement of the differential equation, we chose to specify fn simply as x rather than x[t].  The solution is then returned as a replacement rule for a pure function that can be used with either symbolic or numerical arguments.  Furthermore, it is not necessary for the symbolic argument to be the same as the dummy argument, labelled t, used to represent the independent variable.  Thus, we obtain functions for the displacement and velocity

{x[τ], (x^′)[τ]}  /. 
      sol1[[1]]

as expressions in terms of τ or

{x[1.5/(Subscript[ω, 0])], (x^′)[1.5/(Subscript[ω, 0])]} /. 
      sol1[[1]]

as numerical results quite easily.

Often one sees DSolve used to obtain a solution for x[t] as follows.

sol2  =  
      DSolve[(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x[t],  t]

In this example we specified fn as x[t] and hence obtained a replacement rule for x[t].

x[t]  /. sol2[[1]]

However, x[1.5/(Subscript[ω, 0])] is not recognized as being related to x[t] and hence cannot be evaluated by using the replacement rule directly.

x[1.5/(Subscript[ω, 0])]  /. 
      sol2[[1]]

To evaluate this quantity using the replacement for x[t] would require a couple of steps.

x[t]  /. sol2[[1]]  /. 
      t -> 1.5/(Subscript[ω, 0])

Nor can we obtain the velocity

(x^′)[t] /. 
      sol2[[1]]

with substitution for x[t] before differentiation.

D[x[t] /. sol2[[1]], t]

It is usually easier to express solutions as pure functions rather than in terms of explicit arguments.

Undetermined constants

If the equations do not completely determine the solution, DSolve produces constants of integration indexed sequentially.  By default, the constants are represented by C[1], C[2], ..., C[n], but the function applied to each index can be changed using the option GeneratedParameters. (It was DSolveConstants before version 5.)  By default, GeneratedParameters->C but C remains a reserved symbol even if the option is changed.

sol3 = DSolve[(x^″)[t] ==  -(Subscript[ω, 0]^2) 
              x[t], x, t, GeneratedParameters->a]

C = 1

C[1]  =  0

a[1]  =  0

x[t]  /. sol3

Initial or boundary conditions

Initial conditions or boundary conditions should be specified in the form of equations also.  DSolve[eqns,fn,var] can then be used to solve a list of differential equations, eqns.  If eqns is sufficiently complete, there will be no undetermined constants left.  For example, if we specify that the system is released from initial displacement (Subscript[x, 0]) with zero initial velocity, we obtain a complete solution.

sol4  =  
      DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
          x[0] == a, (x^′)[0] == 0},  x,  
        t]

Ordinarily, we would prefer to specify the initial displacement as (Subscript[x, 0]), but since confusion can arise between x and (Subscript[x, 0]), Mathematica simple echoes the input without producing a solution.

DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x[0] == (Subscript[x, 0]), (x^′)[0] == 0},  x,  
      t]

This problem is solved by using the Symbolize command provided by the Notation package to identify (Subscript[x, 0]) as a symbol that is distinct from x.  We then obtain a solution of the desired form.  (However, you must remember that Symbolize can only be invoked from the Notation palette and will not work properly if typed directly into an input cell.)

Symbolize[(Subscript[x, 0])]

DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x[0] == (Subscript[x, 0]), (x^′)[0] == 0},  x,  
      t]

More generally, if we symbolize the initial velocity also, we can produce a general solution for arbitrary initial conditions.

Symbolize[(Subscript[v, 0])]

shm  = 
      DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
          x[0] == (Subscript[x, 0]), (x^′)[0] == (Subscript[v, 0])},  x,  
        t]

{x[t], (x^′)[t]} /. 
      shm[[1]]

Alternatively, we can present boundary conditions in terms of the solutions for several independent values of the independent variable.  For example, in the case of simple harmonic motion, specification of the displacements at two times is sufficient to determine both constants of integration.

Symbolize[(Subscript[x, 1])]

DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) x[t],  
        x[0] == (Subscript[x, 0]),  x[1/(Subscript[ω, 0])] == (Subscript[x, 1])},  x,  
      t]

Although the initial conditions are usually specified in terms of simple equations like x[0] == (Subscript[x, 0]), it is also possible to use more complicated equations whose solutions may include several sets of  initial conditions, in which case DSolve produces several solutions to the differential equation.  The following variation of the shm problem appears contrived, but illustrates the idea.

DSolve[{(x^″)[t] ==  -(Subscript[ω, 0]^2) 
            x[t],  ((((Subscript[ω, 0]) x[0] - (x^′)[0]))^2) == ((Subscript[v, 0])^2),  (Subscript[ω, 0]) x[0] + (x^′)[0] == 0},  
      x,  t]

Of course, you will probably need some subsidiary condition to choose between the solutions produced.

What can it do?

DSolve claims to be able to solve any linear differential equation with constant coefficients, regardless of order, and many linear equations with nonconstant coefficients up to second order.  Many common special functions were originally defined as solutions to differential equations arising in diverse problems in physics, and hence those equations will usually be recognized.

DSolve[(y^″)[x]  -  x y[x] == 0, y[x], 
      x]

DSolve[(z^2) (y^′)′[z]+z y'[z]+((z^2) - (n^2)) y[z]==0,y,z]

DSolve[ (y^′)′[z] + (1 - (n/(z^2))) y[z] == 
        0, y, z]

However, many relatively simple differential equations cannot be solved in terms of named functions.

DSolve[ (y^′)′[z] + Sin[z] y[z] == 0, y, 
      z]

DSolve[ (y^′)′[z] + (1 - (n/(z^n))) y[z] == 
        0, y, z]

Solve (y^′)  +  y tan(x)  =  sin(2  x) given y = 2 when x = 0 and verify the solution.

Obtain the symbolic solution to a second-order linear differential equation with constant coefficients in the general form a (y^″)  +  b (y^′)  +  
          c y  =  0 and interpret the solution.

Solve ((((y^′)))^2) +  y = a and verify the solution.

Solve (y^′) +  (y^2) = a and verify the solution.

Numerical solution

Basic syntax

The command for numerical solution of differential equations is NDSolve with the basic syntax  NDSolve[eqns,fn,{var,min,max}], where eqns must now be complete enough for a proper numerical solution and var must be given a definite range between min and max.  Hence, eqns must specify a sufficient number of initial conditions and all symbols must evaluate to numbers.  Although the initial conditions usually specify values and derivatives directly, more complicated relationships can be employed provided that the number of conditions is sufficient.  These conditions may be expressed in terms of points outside the range of the solution sought.  If the equations for the initial conditions admit multiple solutions, NDSolve will produce multiple solutions also.

NDSolve returns its results in the form of interpolating functions.  For example, consider the equation of motion for an ordinary pendulum without making the simplifying assumption of small-amplitude oscillations.

ndsol1  = 
      NDSolve[{(θ^″)[
              t] ==  -Sin[θ[t]], θ[
              0] == ((3  π)/4), (θ^′)[0] == 
            0},  θ,  {t, 0, 10}]

An InterpolatingFunction is used to represent an approximate function whose values are obtained by interpolation within a table.  It has the form of a pure function.  Thus, InterpolatingFunction[domain,table][x] computes the function at argument x using an interpolating function defined within a specified domain based upon values contained in a table.  In standard output only the domain is printed, with the table indicated schematically by <>.  Note that if the argument is outside the domain, a warning message is issued and an extrapolated value returned; however, extrapolation tends to go bad very quickly outside the domain of the interpolating function and should be avoided.

θ[t] /. 
      ndsol1[[1]]

An interpolating function can be used like any other function.  In particular, its derivative can be evaluated using D or Derivative.  Thus, it is a simple matter to construct the phase plot for the large-amplitude pendulum using ParametricPlot to plot the angular velocity against the angular displacement.  Instead of the simple circle characteristic of the simple harmonic oscillator, we now obtain a flattened curve that resembles an American football.

ParametricPlot[
        Evaluate[{θ[t], (θ^′)[t]} /. 
            ndsol1[[1]]], {t, 0, 10}, 
        AspectRatio -> Automatic];

Driven pendulum

(((∂^2) θ)/(∂(t^2)))+γ ((∂θ)/(∂t))+ ((Subscript[ω, 0])^2) sin(θ)==f cos(ω t)

Systems of equations

NDSolve will also solve systems of differential equations.  For example, another method for producing a phase plot is to reduce a second-order differential equation to a pair of first-order equations and then to solve for the displacement and velocity as if these were independent functions.

ndsol3  = 
      NDSolve[{ω[t]  ==  (θ^′)[
              t], (ω^′)[
              t] ==  -Sin[θ[t]], θ[
              0] == ((3  π)/4), ω[0] == 
            0},  {θ, ω},  {t, 0, 10}]

Although there is only one solution, that solution consists of two replacement rules.

ParametricPlot[
        Evaluate[{θ[t], ω[t]} /. 
            ndsol3[[1]]], {t, 0, 10}, 
        AspectRatio -> Automatic];

Another interesting system of equations are the Lorenz equations, which arise in the theory of convection.  With appropriate parameters and initial conditions, one obtains the famous butterfly attractor.

\(\(lorenzEqs\  = \ {
		\(x\^′\)[
              t]\  == \ α \((\ \(-x[t]\) + 
                  y[t])\), 
		\(y\^′\)[
              t]\  == \ \(-x[t]\)\ z[t]\  + \ β\ x[t]\  - \ 
              y[t], 
		\(z\^′\)[t]\  == \ 
            x[t]\ y[t]\  - \ γ\ z[t]\ };\)\)

lorenzPlot[values_, ic_, tmax_] :=  
      ParametricPlot3D[
        Evaluate[{x[t], y[t], 
              z[t]} /. NDSolve[
                Join[lorenzEqs /. values, ic], {x, y, z}, {t, 0, tmax}, 
                MaxSteps -> 
                  10000][[1]]], {t, 0, 
          tmax}, PlotPoints -> 3000]

A trajectory using the customary parameters illustrates the complexity of this motion.

lorenzPlot[{α -> 10, β -> 28, γ -> 
            (8/3)}, {x[0] == 0, y[0] == 1, z[0] == -1}, 
        50];

Sensitivity to initial conditions

Accuracy issues

NDSolve employs several sophisticated adaptive algorithms which attempt to optimize both speed and accuracy without user intervention, but no algorithm can succeed for every conceivable set of differential equations.  Consequently, a vast literature has accumulated over the last 100 years or so discussing a large variety of numerical methods and comparing their strengths, weaknesses, and realms of applicability.  Although it is certainly beyond the scope of this course to review the literature on the subject or to implement many of these competing algorithms, one must be aware of the options available for regulating the accuracy of the standard algorithm, for Murphy's law decrees that your first interesting problem will challenge the default algorithm.  At that time you may be motivated to produce a custom integrator suitable to your application or may seek an alternative package.  In either case, one must have some insight into the nature of the equation at hand and the special challenges it might present.  Nevertheless, after all of these obligatory caveats, one should recognize that the default algorithm is already as good or better than most of the alternatives and that if your problem looks intractable in its present form the best approach is often to seek a better representation of the problem than a better numerical algorithm — brute computation is a poor substitute for insight.

A very good discussion of the accuracy options available to NDSolve can be found in Section 3.9.7 of the Mathematica book.  First, we list those options.

Options[NDSolve]

AccuracyGoal specifies the number of decimal digits for the solution itself, while PrecisionGoal governs the relative precision of the solution.  If the solution varies strongly in magnitude, the relative precision may be most important, particularly if small errors in small values tend to grow into large errors elsewhere, which is a property of unstable systems.  Setting AccuracyGoal->Infinity instructs the program to test PrecisionGoal only.  With automatic settings, AccuracyGoal and PrecisionGoal are about half the WorkingPrecision; hence, increasing either goal generally requires increasing WorkingPrecision also.  If necessary, NDSolve adjusts its step sizes to meet these goals provided that MaxSteps is not exceeded.  MaxSteps may be exceeded either when there are accuracy problems or when the range of integration is very large, in which case one can increase it, perhaps to Infinity if you are patient and brave.  Occasionally it is necessary to reduce StartingStepSize if the solution varies rapidly near its starting point.

Consider, for example, the following simple equation.

(y^″)[x]  -  (λ^2) 
          y[x]  ==  0,  
    y[0] == 
      1,  (y^′)[
          0] == -λ   ⟹  y[
              x]  =  (e^(-λ x))

Although the symbolic solution is trivial, numerical algorithms may encounter considerable difficulty if λ happens to be large.  The reason for this behavior is easy to explain for such a simple equation.  Consider the original equation without the initial conditions.

(y^″)[x]  -  (λ^2) 
            y[x]  ==  0   ⟹  y[x]  =  
      A (e^(-λ x))  +  
        B (e^(λ x))

The second-order differential equation has two linearly independent solutions whose mixing is determined by the initial conditions.  Although we chose initial conditions which selected the exponentially decaying solution and eliminated the growing solution, numerical methods find it much more difficult to suppress the divergent solution.  Any numerical error in the initial conditions or the integration procedure will introduce a small admixture of the spurious solution, which then grows uncontrollably until it overwhelms the correct solution.  The present equation is described as stiff because the coefficient that multiplies errors is large, here (λ^2), so that small errors are quickly amplified.

The following functions can be used to compare numerical solution with suitable acccuracy option with the true solution.

stiffEq[λ_] := {(y^″)[
                x] - (λ^2) y[x] == 0, 
          y[0] == 
            1, (y^′)[
              0] == -λ};
stiffSol[λ_, xmax_, opts___Rule] := 
        y[x] /. NDSolve[stiffEq[λ], y, {x, 0, xmax}, 
              opts][[1]];
trueSol[λ_] := 
      y[x] /. DSolve[stiffEq[λ], y, x]

Consider the two plots show below.  Both solve the same equation with the same large λ, but the first uses a larger range of x than the second even though  both are plotted for the same range.  The first solution appears to start plausibly, but a spurious divergence occurs for large x that is absent in the second solution with a more restricted range.  Evidently, a larger step size is used for the larger range and numerical errors near the starting point produce an intolerably large admixture of the spurious divergent solution.

Plot[Evaluate[stiffSol[20, 30]], {x, 0, 10}, 
        PlotRange -> {Automatic, {0, 1}}];
Plot[Evaluate[stiffSol[20, 10]], {x, 0, 10}, 
        PlotRange -> {Automatic, {0, 1}}];

A small improvement in the range of validity can be obtained by reducing the initial step size

DisplayTogether[
        LogPlot[Evaluate[
            Abs[stiffSol[20, 10, StartingStepSize -> (10^(-6))]]], {x, 
            0, 2}, PlotRange -> {Automatic, {(10^(-9)), (10^1)}}], 
        LogPlot[Evaluate[trueSol[20]], {x, 0, 2}, 
          PlotStyle -> {{RGBColor[1, 0, 0], 
                Dashing[{0.02, 0.02}]}}]];

or by increasing the accuracy goal

DisplayTogether[
        LogPlot[Evaluate[
            Abs[stiffSol[20, 10, AccuracyGoal -> 10, 
                WorkingPrecision -> 20]]], {x, 0, 10}, 
          PlotRange -> {Automatic, {(10^(-30)), (10^1)}}], 
        LogPlot[Evaluate[trueSol[20]], {x, 0, 10}, 
          PlotStyle -> {{RGBColor[1, 0, 0], 
                Dashing[{0.02, 0.02}]}}]];

or by increasing the precision goal

DisplayTogether[
LogPlot[
          Evaluate[
            Abs[stiffSol[20, 10, PrecisionGoal -> 10, 
                WorkingPrecision -> 20]]], {x, 0, 10}, 
          PlotRange -> {Automatic, {(10^(-25)), 
                (10^1)}}], 
LogPlot[
          Evaluate[
            Abs[stiffSol[20, 10, PrecisionGoal -> 15, 
                WorkingPrecision -> 30]]], {x, 0, 10}, 
          PlotRange -> {Automatic, {(10^(-25)), (10^1)}}, 
          PlotStyle -> {{RGBColor[0, 0, 
                  1]}}], 
LogPlot[
          Evaluate[trueSol[20]], {x, 0, 10}, 
          PlotStyle -> {{RGBColor[1, 0, 0], 
                Dashing[{0.02, 0.02}]}}]];

Notice that the working precision must be about twice the accuracy or precision goal for these to be completely effective.  When a symbolic solution is not available, one can compare results for two values of the accuracy or precision goals to evaluate the probable range of validity for a numerical solution.  However, if one really needs to push the limits of numerical methods, it is probably better to reformulate the problem!

DisplayTogether[
LogPlot[
          Evaluate[
            Abs[1 - stiffSol[20, 10, PrecisionGoal -> 10, 
                    WorkingPrecision -> 20]/
                  stiffSol[20, 10, PrecisionGoal -> 20, 
                    WorkingPrecision -> 40]]], {x, 0, 2}, 
          PlotRange -> {Automatic, {(10^(-10)), 
                (10^1)}}], 
LogPlot[
          Evaluate[
            Abs[1 - stiffSol[20, 10]/
                  stiffSol[20, 10, PrecisionGoal -> 20, 
                    WorkingPrecision -> 40]]], {x, 0, 2}, 
          PlotRange -> {Automatic, {(10^(-10)), (10^1)}}]];

Although this is a trivial example, nontrivial examples of a less contrived nature often present similar problems and require more detailed analysis to isolate their causes.  Fortunately, the accuracy options provided by NDSolve will often suffice where other standard algorithms fail, but it doesn't pay to become complacent — it is usually worthwhile to investigate more carefully and to gain more insight into the behavior of the system.  Sometimes rapid change is an important feature of the system and other times it may be an artifact of numerical instability, but it is your responsibility to distinguish between these possibilities!

A stiff equation

L (((∂^2) q)/(∂(t^2)))+R ((∂q)/(∂t))+(q/C)==(Subscript[V, 0]) cos(ω t)

(((∂^2) x)/(∂(t^2)))+γ ((∂x)/(∂t))+ ((Subscript[ω, 0])^2) x (1 - α (x^2)) ==f cos(ω t)

-(ψ^″)[x]  +  ν[
            x] ψ[x]  ==  ψ[x]

(y^′)[x]  +  x y[x]^2  ==  
      a,   y[0] == 0