Symbolic Differential Equations

Advances in DSolve

Devendra Kapadia,

Wolfram Research, Inc.

In this talk, I would like to do the following :

Give a quick tour of  DSolve

Discuss the Cauchy Problem for First-Order PDEs

Show some solutions of ODEs using the new special functions

Speak about dummy variables in DSolve and other improvements

Highlight the role of the Painleve Equations in DSolve

At present, DSolve can return exact solutions for :

-  Ordinary differential equations (linear, nonlinear, single, system)

-  Partial  differential equations  (general first-order PDEs, some higher order equations)

-  Differential - algebraic equations (linear constant coefficient systems)

This ODE is a model for an aging spring. It has an exponential function in its coefficient and is solved by transforming it to Bessel's equation.

AgingSpring = x''[t] == - 4*Exp[-t/4]*x[t];

sol=DSolve[{AgingSpring , x[0]==1, x'[0]==2}, x,t]

Plot[x[t]/.sol[[1]], {t,0,30}, PlotStyle -> {Thickness[0.0075],Red}, PlotPoints -> 100, ImageSize -> 400]

Consider the following linear first-order PDE.

LinearPDE = y*D[u[x, y], x] + x*D[u[x, y], y] == -4*x*y*u[x, y];

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

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

We obtain a Gaussian function by making a particular choice of the function C[1].

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

The solution surface may be plotted as follows :

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

We begin by setting up the problem using a differential equation, an algebraic constraint and an initial condition.

DifferentialEquation = x'[t] - y[t]==Cos[t];

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

InitialCondition = x[Pi/2]==3;

DAE= {DifferentialEquation, AlgebraicConstraint, InitialCondition};

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

Plot[{Tooltip[x[t]/.sol, Style["x[t]", Red,18]],Tooltip[ y[t]/.sol, Style["y[t]", Blue,18]],Tooltip[ (x[t]+ y[t])/.sol, Style["x[t]+ y[t]", RGBColor[0,0.7, 0.4],18]]}, {t,3,12}, PlotStyle -> {{Thickness[0.01],Red}, {Thickness[0.01],Blue}, {Thickness[0.01],Dashing[{0.05}],Green}}]

A first - order partial differential equation with unknown function u(x,y) is :

1. Linear if it it is linear in  u, ((∂u)/(∂x)) and ((∂u)/(∂y)). For example,

TransportEquation =((∂u(x, y))/(∂x))+ c ((∂u(x,y))/(∂y))==0

2. Quasilinear if it it is linear only in ((∂u)/(∂x)) and ((∂u)/(∂y)). A good example is

BurgersEquation =   ((∂u(x,y))/(∂x))+ u(x,y)((∂u(x,y))/(∂y))==0

3. Nonlinear if it is not of type 1 or 2 above. For example,

EikonalEquation =    (((((∂u(x,y))/(∂x))))^2)+ (((((∂u(x,y))/(∂y))))^2)==1

The nature of the PDE (linear, quasilinear, nonlinear) determines the type of solutions returned by DSolve.

In the  linear case, we return the general solution of the PDE.

SetOptions[DSolve, GeneratedParameters ->( Subscript[f, #]&)];

SolutionForTransportEquation = 
  DSolve[D[u[x, y], x] + c*D[u[x, y], y] == 0, u[x, y], {x, y}]

For  quasilinear equations, we again return the general solution, but it may be in implicit form.

SolutionForBurgersEquation = 
  DSolve[D[u[x, y], x] + D[u[x, y], y]*u[x, y] == 0, u[x, y], {x, y}]

Nonlinear first-order PDEs are solved in terms of complete integrals, depending on arbitrary constants.

SetOptions[DSolve, GeneratedParameters ->( Subscript[c, #]&)];

EikonalSolution = DSolve[
   D[u[x, y], x]^2 + D[u[x, y], y]^2 == 1, u[x, y], {x, y}]

Given  values for u[x, y]  on a curve  in the x - y plane (say, u[x, 0] == f[x]), find a solution of the PDE which also satisfies these initial conditions.

The problem is usually solved by reduction to a system of ODEs along the so - called characteristic curves.

Plot[{Tooltip[1/t^2,Style[ "Initial Curve", {Red ,Bold,18}]], Tooltip[-3+7 t, Style["Characteristic Curve", {Blue ,Bold,18}]],Tooltip[-2+7 t,Style[ "Another Characteristic Curve", {Blue ,Bold,18}]],Tooltip[-1+7 t,Style[ "Yet Another Characteristic Curve", {Blue ,Bold,18}]],7t, 1+7 t}, {t,0,1}, Ticks -> None,PlotStyle ->{{Thickness[0.007], Red}, Blue, Blue, Blue, Blue, Blue, Blue, Blue}]

In  Mathematica 6.0, we either use the general solution or the complete integral or the envelopes of complete integrals to solve the Cauchy problem.

Here is the solution of the transport equation for arbitrary initial conditions.

TransportEquation = D[u[x, y], x] +  D[u[x, y], y] == 0;

InitialCondition = u[x, 0] == f[x];

CauchyProblem = {TransportEquation, 
    InitialCondition};

DSolve[CauchyProblem, u[x, y], {x, y}]

The general solution for Burger' s Equation is in implicit form.

However,  the initial - value problem for this equation can sometimes be solved explicitly.

BurgersEquation = D[u[x, y], x] +  u[x, y]*D[u[x, y], y] == 0;

InitialCondition = u[x, 0] == 1/(x + 1);

DSolve[{BurgersEquation, InitialCondition}, u, {x, y}]

The following solves a Cauchy problem for the Clairaut equation.

ClairautEquation = u[x, y] == x*D[u[x, y], x] + y*D[u[x, y], y] + (1/2)*(D[u[x, y], x]^2 + 
       D[u[x, y], y]^2);

InitialCondition = u[x, -(x^2/2)] == 1/2;

DSolve[{ClairautEquation ,InitialCondition}, u, {x, y}]

{ClairautEquation ,InitialCondition}/.%[[1]]//Simplify

In a recent paper, Herbert Wilf showed how to find the generating functions of the Binomial and Stirling numbers by solving a Cauchy problem.

In V 6.0, we can use DSolve for this purpose as shown in the next example.

CauchyProblemForStirlingS2 = 
   {D[u[x, y], x] == y*D[u[x, y], y] + y*u[x, y], u[0, y] == 1};

The generating function for StirlingS2[n, k] can be now be obtained as follows:

GeneratorForStirlingS2[x_, y_] =Simplify[ u[x,y]/.DSolve[CauchyProblemForStirlingS2, u, {x,y}][[1]]]

The value of StirlingS2[11, 5] (say) can now be calculated as follows.

SeriesCoefficient[Series[GeneratorForStirlingS2[x,y],{x,0,11},{y,0,5}],{11, 5}]*11!

StirlingS2[11,5]

The details are available in a forthcoming paper, "The Cauchy Problem for First-Order Partial Differential Equations", (to be published in The Mathematica Journal).

ODEs for the New Special Functions

There are several new special functions in Mathematica 6.0 - ParabolicCylinderD, SpheroidalPS, KelvinBei, etc.

Plot[{Tooltip[ParabolicCylinderD[3., x],Style[ "ParabolicCylinderD", {Red ,Bold,18}]], Tooltip[KelvinBei[2.,x],Style[ "KelvinBei", {Blue ,Bold,18}]], Tooltip[SpheroidalPS[2,0,1/3,x],Style[ "SpheroidalPS", {Brown ,Bold,18}]]}, {x,-1,4}, PlotStyle -> {{Thickness[0.01], Red}, {Thickness[0.01], Blue}, {Thickness[0.01], Brown}}]

For this example, we get the solution given in  Whittaker and Watson (pg 347, Section 16.51).

eqn =  y''[x] + ((1/2) + n - x^2/4)*y[x];

DSolve[eqn== 0, y, x]

eqn/.%[[1]]/.{x-> Random[], n-> 11/3}//Simplify//Chop

This solution for the basic spheroidal equation is new in V 6.0.

eqn1=((g^2)-(g^2) (x^2)-((m^2)/(1-(x^2)))+SpheroidalEigenvalue[n,m,g]) y[x]-2 x y'[x]+(1-(x^2)) y''[x];

sol = DSolve[eqn1 == 0, y, x]

FullSimplify[eqn1 /. sol[[1]]]

The four Kelvin functions satisfy a fourth - orde linear ODE.

eqn =(a* x+b)^4*(1/a^4)*y''''[x]+2(a*x+b)^3*(1/a^3)*y'''[x]-(1+2 v^2)((a*x+b)^2*(1/a^2)*y''[x]-(a*x+b)*(1/a)*y'[x])+(v^4-4v^2+(a*x+b)^4)y[x];

DSolve[eqn==0, y,x]

eqn/.%[[1]]/. {x-> N[(37/20), 20], v-> 4, a-> 5, b-> -2}//Simplify//Chop

The dummy variables in unevaluated integrals are now labelled K[1], K[2] etc. This replaces module K$ variables in earlier versions.

DSolve[y'[x] == f[x], y, x]

DSolve[y''[x] == f[x], y, x]

In V 6.0, we can recognize a larger number of basic equations such as the Struve equation below, so that the most appropriate solution can be returned.

Clear[y,z]

DSolve[z^2 y''[z]+z y'[z]+(z^2-n^2) y[z]==2/Pi z^(n+1)/(2 n-1)!!,y,z]

Here is a travelling wave solution for the Korteweg-deVries  equation.

SetOptions[DSolve, GeneratedParameters -> C];

KdV = {D[u[x, t], t] + D[u[x, t], {x, 3}] + 6*u[x, t]*D[u[x, t], x] == 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 -> Orange, FillingStyle -> GrayLevel[0.6], Filling -> Axis, ImageSize-> 800], {a,0.1,3,0.02}, DefaultDuration -> 30]

The solutions of ODES may have fixed or movable singularities as seen in the following examples.

DSolve[(x-3) y'[x]+y[x]==0,y[x],x]   (* fixed pole *)

DSolve[y'[x] +y[x]^2 ==0, y[x], x]   (* movable pole *)

DSolve[y'[x] +y[x]^3==0, y[x], x][[1]]  (* movable algebraic branch point *)

A critical point of a differential equation is one where the solution can become multiple-valued.

Around 1900, the French mathematician Paul Painlevé showed that  ODEs of type

y''(x)==R(x,y(x),y'(x))

with no movable critical points, where R is  a rational function  of its arguments, can  be reduced to fifty canonical equations.

These fifty equations are called the Painlevé equations and are divided into two groups :

1. Forty - four of these equations can either be solved explicitly (soluble) or have first integrals (semi - transcendental).

2. The remaining six equations require new special functions and are called Painlevé transcendents.

The  Painlevé transcendents are labelled as Painlevé  I - VI.

Painlevé  Equation 3 (solution is given in terms of Weierstrassian elliptic functions):

DSolve[y''[x]==6 ((y[x])^2)-(1/24),y[x],x]

Painlevé  Equation 17 (solution is given in terms of polynomials for positive integers m):

DSolve[w''[z]==(((m-1) ((w'[z])^2))/(m w[z])),w,z]

m=5;

DSolve[w''[z]==(((m-1) ((w'[z])^2))/(m w[z])),w,z]

Clear[m]

The second Painleve transcendent is given by

PainleveII=y''[x]==2 ((y[x])^3)+x y[x]+α;

We can find special rational solutions of this equation using the family of Vorob'ev - Yablonskii polynomials Q[n, z], given by

Clear[Q]

Q[0,x_]:=1

Q[1,x_]:=x

Q[n_,x_]:=Q[n,x]=Together[((x ((Q[n-1,x])^2)-4 Q[n-1,x] (∂_({x,2}))Q[n-1,x]+4 ((((∂_x)Q[n-1,x]))^2))/(Q[n-2,x]))]

Q[6, x]

The roots of these polynomials lie on approximate equilateral triangles.

sol =NRoots[Together[Q[14, x]]==0, x];ListPlot[{Re[x], Im[x]}/.{ToRules[sol]}, PlotStyle -> {PointSize[0.03], Red}, Background -> Cyan, AspectRatio -> 1]

The rational solutions of Painlevé II can be constructed using the Vorob' ev - Yablonskii polynomials as follows :

y2[x_, 0]:=0

y2[x_, α_]:= D[Log[Q[α-1, x]/Q[α, x]], x]

Table[{α, y2[x,α]}, {α, 1, 4}]//Apart//TableForm

An elementary exposition of the Painlevé equations using Mathematica 6.0  is  available in notebook form
(send email to dkapadia@wolfram.com).

The following new features of DSolve make it more powerful and easier to use in V 6.0.

Solution of initial - value problems for first - order PDEs.

Support for the differential equations associated with the new special functions.

Other improvements such as the simpler dummy variables in unsolved integrals.

In future versions, we plan to do the following :

Implement solutions for the boundary value problems associated with the classical second - order PDEs.

Develop algorithms for new classes of ODEs  such as the Lamé and Heun equations.