Symbolic Differential Equation Solving

Advances in DSolve

Devendra Kapadia

Wolfram Research, Inc.

Show how to perform some basic tasks using DSolve.

Take a trip through the Advanced Documentation.

Discuss  factorization methods for differential operators.

Demonstrate a travelling wave solution for the KdV equation.

Give some insight into the working of DSolve.

The following equation is an example of Mathieu's equation.

Mathieu = y''[x] + (2 Cos[2x] + 7)y[x] == 0;

Here is the general solution of this equation.

sol =DSolve[Mathieu, y, x]

Mathieu/.sol//Expand

We can plot the solution for specific values of C[1] and C[2].

Plot[y[x]/.sol/.{C[1]-> 1, C[2]-> 3/4}, {x, -5,5}, PlotStyle ->{Thickness[0.01],  Orange}]

The following equations describe a plane curve through the origin whose curvature is equal to the parameter value 's' at 
every point.

PlaneCurve = {x'[s]==Cos[t[s]], y'[s]==Sin[t[s]], t'[s] ==s, x[0]==0, y[0]==0, t[0]==0};

The solution of the ODEs contains Fresnel functions.

sol =DSolve[PlaneCurve, {x,y, t}, s]

A graph of the curve clearly shows that the curvature becomes large as 's' approaches infinity.

ParametricPlot[{x[s], y[s]}/.sol[[1]], {s,-30,30}, PlotStyle -> Red]

The following is  a simple first order linear PDE.

z:=u[x,y]

p:=D[u[x,y],x]

q:=D[u[x,y],y]

LinearPDE=y*p+x*q ==-4x*y*z;

The general solution of this PDE contains an arbitrary function C[1].

sol=DSolve[LinearPDE, u, {x,y}]

This checks that the solution is correct.

LinearPDE/.sol[[1]]//Simplify

Finally, we obtain a Gaussian function by making a particular choice of C[1].

particularsolution = (u[x,y]/.sol[[1]]/.C[1][a_]:>ExpandAll[Exp[-2a]])

Plot3D[particularsolution, {x,-2, 2}, {y,-2,2}, PlotPoints -> 30,ColorFunction -> Hue, PlotRange -> All]

In this example we use DSolve to obtain a function whose gradient vector is known.

GradientVector ={p==2x,q==2y};

DSolve[GradientVector, z, {x,y}]

Clear[z,p,q]

Differential = x'[t] - y[t]==Sin[t];

Algebraic  = x[t]+ y[t] ==1;

Initial = x[Pi/2]==1;

DAE= {Differential, Algebraic, Initial};

sol=DSolve[DAE, {x,y}, t]

Plot[{x[t]/.sol, y[t]/.sol,(x[t]+ y[t])/.sol}, {t,3,12}, PlotStyle -> {{Thickness[0.01],Red}, {Thickness[0.01],Blue}, {Thickness[0.01],Dashing[{0.05}],Green}}]

Goals of the Advanced Documentation

- Provide a catalog for the equations which can be handled by DSolve.

- Give a brief description of the algorithms used for solving various classes of equations.

- Include a self-contained introduction for beginning users with tips for applying DSolve effectively.

Some salient features

- There is a table with hyperlinks and historical data for many types of ODEs.

- The classification of PDEs is discussed in detail.

- Examples are provided for special kinds of input (such as BVPs with Piecewise coefficients).

Working with DSolve

- Designed to be a comprehensive introduction to DSolve functionality.

- Discusses questions such as the use of inexact numbers in DSolve.

- Emphasizes the need for checking that your problem is well-posed.

The notion of a symmetric power

Suppose  {u, v} is a basis for the solution space  of a homogeneous linear second order ODE.

For example, if we use Airy's equation,  then u = AiryAi[x] and v = AiryBi[x].

Airy= y''[x]-x*y[x];

DSolve[Airy==0, y[x], x]

Then the second symmetric power of this equation is a linear third order ODE with  solution basis given by {(u^2),u v,(v^2)}.

More generally, the nth symmetric power of the second order ODE is a linear (n+1)th order equation with solution basis given by {(u^n), (u^(n-1)) v,..., (v^n)}.

The Bronstein-Mulders-Weil-van Hoeij algorithm

1.  Mark van Hoeij provided an easy way of calculating the coefficients of the underlying second order ODE using
      only the top two coefficients of the higher order ODE. This assumes that the input is a symmetric power.

2.  If van Hoeij's criterion is satisfied,  the  (n-1)th symmetric power of the second order ODE is calculated using an 
     iterative method due to Bronstein, Mulders and Weil.

Here, the internal (undocumented) symbol DSolve`DSolveSymmetricPower is used to calculate the eigth symmetric power of the Airy equation.

Clear[x, y,z]

DSolve`DSolveSymmetricPower[0, -x, 9, y,x]//Expand

3. We calculate the general solution of the second order ODE and use it to write down the general solution of the higher 
    order  ODE.

An example of order 10

Here we solve a tenth order ODE in terms of Airy functions.

Airy10 = (DSolve`DSolveSymmetricPower[0, -x, 10, y,x]//Expand)

(sol =DSolve[Airy10==0, y,x])//Timing

Airy10/.sol[[1]]//Expand

Basic definition of a symmetric product

Suppose that we are given two second order ODEs with solution bases {f, g} and {r, s}.

A fourth order ODE is a symmetric product of these second order ODEs if  {f r, f s, g r, g s} is a basis 
          for its general solution.

The key technical tool  used in  van Hoeij's algorithm

The algorithm for solving symmetric products is due to van Hoeij and relies on finding rational-exponential 
          solutions of an associated third order ODE.

In the following example, the (undocumented) symbol DSolve`DSolveExponentialSolutions picks out the 
           exponential solution (e^(2 x)).

DSolve[(-2+2x) y[x]+(-1-x) y'[x]+y''[x]==0, y, x]

DSolve`DSolveExponentialSolutions[{-2+2x, -1-x, 1}, 2,x]

Since this functionality has been available since Mathematica 4.1, the implementation of van Hoeij's 
           algorithm was relatively straightforward.

An example involving Airy and Hermite functions

The following fourth order ODE is a symmetric product of the equations for the Hermite and Airy functions.

DSolve[(((4+5 x+4 (x^2)+7 (x^3)+6 (x^4)+3 (x^5)+(x^6)) y[x])/(1+x+(x^2)))+(((-5+x+5 (x^2)+2 (x^3)) y'[x])/(1+x+(x^2)))+(2-2 x+2 (x^2)) y''[x]+(((-1-2 x) Derivative[3][y][x])/(1+x+(x^2)))+Derivative[4][y][x]==0, y,x]//Timing

Building an ODE from lower order factors

Starting with two lower order  differential operators, it is possible to build a higher order operator by
       composition as shown in the following example.

LeftFactor=(D[#,x]-#)&; (* y'[x]-y[x]*)

RightFactor=y''[x]-x*y[x];(* Airy's equation *)

ThirdOrderOperator=Collect[LeftFactor[RightFactor],{y[x],y'[x],y''[x],y'''[x]}]

The problem of factorization is to find the lower order factors, if any, of a given higher order operator.

Outline of the algorithm

The algorithm relies on modern versions of the associated equations method, due to F. Schwarz and M. Bronstein.

Step 1: Find the left factor L and the right factor R, using the  (undocumented) symbol   
                 DSolve`DSolveFactorLinearDifferentialOperator.

DSolve`DSolveFactorLinearDifferentialOperator[{-1+x, -x, -1, 1}, 3,y,x]//Timing

Step 2: solve the lower order ODE L = 0.

Step 3: Solve the inhomogeneous ODE R = f[i] , where f[i] are elements of the solution basis 
                  of L   (time-consuming step).

An example due to Bronstein

The following fourth order ODE appears in  the paper by Bronstein and is solved using factorization in Mathematica 5.2.

eqn = y''''[t] -2 t y''[t] -2 y'[t] +t^2 y[t];

(sol=DSolve[eqn==0, y,t])//Timing//InputForm

Here is a plot of the solution for a particular choice of arbitrary constants.

ListLinePlot[Table[{t,Evaluate[y[t]/.sol[[1]]/.{C[1]-> 1, C[2]-> 1/3, C[3]-> 3/4, C[4]-> 1/7}]}, {t, -8.,1.5, 0.1}], PlotStyle -> {Thickness[0.01], Red}]

During the last few decades, there has been great interest (and success) in finding special solutions such as solitons, for nonlinear partial differential equations.

These special solutions belong to a broad class called travelling wave solutions.

If x denotes the space variable and t denotes the time variable, then a travelling wave solution usually occurs as polynomial in a 
hyperbolic function (Tanh or Sech)  (for example, (sech^2)(x (Subscript[c, 1])+t (Subscript[c, 2])+(Subscript[c, 3]))).

We have implemented an algorithm due to Baldwin, Goktas, Hereman and others for finding travelling wave solutions of the Korteweg-deVries (KdV) and other nonlinear PDEs.

This is a very general algorithm and can be used to find solutions in terms of Tanh, Sech, JacobiSN and JacobiCN. It works 
for both single as well as systems of PDEs. However, we will refer to it simply as the "Tanh" method.

The algorithm relies on the fact that the derivatives of these functions are polynomial expressions of the functions. 
For example,

T1=Tanh[x];

D[T1,x]

1-T1^2

Simplify[%]

Clear[T1]

If we look for solutions which are polynomials in T1, then the problem is eventually reduced to solving a nonlinear system of algebraic equations.

Here is a travelling wave solution for the KdV equation.

KdV= {D[u[x, t], t] + 6*u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 3}] == 0};

sol = DSolve[KdV, u, {x, t}]

We can get a (Sech^2) soliton by choosing the parameters carefully, as in the animation below.

Animate[Plot[Evaluate[u[x,a]/.sol[[1]]/.{C[1]-> 1, C[2]-> -4, C[3]-> 1/4}], {x,-7,20}, PlotStyle -> {Thickness[0.01],Blue}, PlotRange -> All, Frame -> True, Background -> Cyan, FillingStyle -> GrayLevel[0.6], Filling -> Axis, ImageSize-> 800], {a,0.1,3,0.02}, DefaultDuration -> 30]

DSolve is equipped with a wide variety of algorithms for solving  differential equations and boundary value problems.

The Advanced Documentation is a useful resource for both beginning and advanced users.

Factorization methods in Mathematica 5.2 add a considerable amount of power to the solution of higher order ODEs.

Travelling wave solutions in closed form are a nice tool for exploring nonlinear partial differential equations.

For future versions, we would like to do the following:

Continue to strengthen the basic ODE solver by adding methods for finding hypergeometric 
and other special function solutions.

More support for linear/nonlinear inital and boundary value problems (first integrals, 
Sturm-Liouville problems).

Add more functionality for solving  PDEs (such as the Cauchy problem).