Symbolic Differential Equations

Advances in DSolve (Mathematica 5.1)

Devendra Kapadia

Wolfram Research, Inc.

Outline of  Presentation

Overview of DSolve (pre - V 5.1)

Abel ODEs with non-constant invariants

Improvements for solution of Hypergeometric type ODEs

Transformations for ODEs with non-rational coefficients

Differential equations with Piecewise coefficients

Complete integrals for first order nonlinear PDEs

Overview of DSolve

The Mathematica function DSolve finds symbolic solutions for the following types of equations.
                            
                            1. Ordinary differential equations (ODEs). 
                            
                            2. Partial differential equations (PDEs).
                            
                             3. Differential-algebraic equations (DAEs).

Here is the solution of Bessel's equation.

bessel  =  
        x^2*y''[x] + x*y'[x]  + (x^2 - 1)*y[x]  == 
          0;

sol = DSolve[bessel,  y, x]

Plot[y[x] /. sol[[1]] /. {C[1] ->  1,  
            C[2] ->  0},  {x,  -15, 15},  
        PlotStyle  ->  {Thickness[0.006], Red}];

Here is the solution for a first order linear PDE.

eqn = 
        x^2*D[u[x, y], x] + 
          y^2*D[u[x, y], y] - (x + y)*u[x, y];

sol = DSolve[eqn == 0,  u,  {x, y}]

Plot3D[u[x, y] /. 
            sol[[1]] /. {C[1][t_] ->  
              Sin[t^2]  + (t/10)},  {x, -5, 5},  {y, -5, 5}, 
        AmbientLight -> Hue[0.7],  
        ImageSize  ->  400];

Here is the solution for a linear DAE with constant coefficients.

daes = {f''[x]  ==  g[x],  f[x]  +  g[x]  ==   x,  
              f[Pi]  ==  1,  f'[Pi]  ==  0};

sol  =  DSolve[daes,  {f,  g},  x]

Plot[{f[x] /. sol,  g[x] /. sol,  
          f[x] + g[x] /. sol},  {x, -5, 5},  
        PlotStyle  -> {{Thickness[0.008], Red},  {Thickness[0.008], 
              Blue}, {Thickness[0.008], Dashing[{0.04}], Hue[0.8]}},  
        ImageSize  ->  400];

Introduction to Abel ODEs

A first order linear ODE has the general form:

y'[x]==P[x] y[x] + Q[x]

Example

DSolve[y'[x]  + y[x]  == 1/x,  y, x]

Note:

1. The general solution of a first order linear ODE will often contain  an unevaluated integral such as 	ExpIntegralEi.

2.  Some first order ODEs are inverse linear (linear after interchanging the roles of x and y).

inverselinearode  =  y'[x] ==  
          1/(-x*y[x]  + Exp[3*y[x]]);

Off[InverseFunction::ifun,  Solve::tdep]

DSolve[inverselinearode,  y[x], x]

A Riccati equation has the form

y'[x]==P[x] y[x] + Q[x] y[x]^2 + R[x]

Example

The solution of the following Riccati equation can be found by using a logarithmic derivative of Airy's equation.

DSolve[y'[x]  +  y[x]^2 == x + 1,  y, x]

An  Abel ODE  is an equation  of the form

y'[x]==f0 + f1 y[x] + f2 y[x]^2 + 
          f3 y[x]^3

where f0, f1, f2 and f3 are functions of x.

Abel ODEs with constant invariant (defined below) can be solved in V 5.0 and earlier.

f0 = 0; f1 = (1/x); f2 = -3; f3 = x;

Invariant = 
      f0 (f3^2) + 
        (1/3) (((2 (f2^3))/9) - f1 f3 f2 - (∂_x) f3 f2 + 
              f3 (∂_x) f2)

AbelODE=y'[x]==f0 + f1 y[x] + f2 y[x]^2 + 
          f3 y[x]^3

sol = DSolve[AbelODE,  y, x]

AbelODE /. sol // Simplify

Clear[f0, f1, f2, f3]

Abel ODEs with non-constant invariants (V 5.1)

In general, it is difficult to solve Abel ODEs with non-constant invariants.

However, if we can find the solution of a particular Abel ODE, then this can be used to generate solutions of all Abel ODEs which are equivalent to it by a coordinate transformation. Thus, the entire equivalence class becomes integrable.

In V 5.1, DSolve can find solutions for many equivalence classes of Abel ODEs with non-constant invariant. The integrable examples are basically of two types.

The ExpIntegralEi term in the solution to this equation comes from solving the linear ODE.

Off[InverseFunction::ifun,  Solve::tdep]

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

The Airy functions in the solution of this Abel ODE  come from the corresponding Riccati equation.

AbelODE=y'[x]==-2 x y[x]^2 + y[x]^3

sol = DSolve[AbelODE, y[x], x]

There is a close connection between solvable classes of Abel ODEs and second order equations of hypergeometric type.

The solutions of this Abel ODE of the second kind are given in terms of Hypergeometic2F1.

AbelODEOfTheSecondKind=y'[x]==((-((2  x)/9) +  x^3 + y[x])/y)[
              x]

sol = DSolve[AbelODEOfTheSecondKind, y[x], x] // Timing

Improvements for solution of Hypergeometric type ODEs

The hypergeoemetric functions play a unifying role in Mathematical Analysis because several elementary functions arise as special cases of them.

Hypergeometric2F1[a, b, c, x] /. {a ->  1/2,  b ->  1,  
        c ->  3/2}

The ODEs for Hypergeometric0F1, Hypergeometric1F1 and Hypergeometric2F1 occur in many applications.

Here is the ODE for Hypergeometric2F1.

Clear[a, b, c, y, x]

sol = 
      DSolve[x (x - 1)*y''[x] +  ((a + b + 1)*x - c)*y'[
                x] + b*a*y[x] == 0,  y[x],  x]

Hence the solution of ODEs of hypergeometric type is of great importance in symbolic computation.

In Mathematica 5.1, several second order linear equations can be solved by transforming them to one of the basic hypergeometric ODEs by a composition of:
                        a. Linear transformations in y[x].   ( y[x] -> P[x] y[x] )
                        b. Bilinear transformations in x.     ( \(x -> \(a\ x\  + \ b\)\/\(\(c\ x\)\(\ \)\(+\)\(\ \)\(d\)\(\ \)\)\) )
                        c. Power transformations in x.       ( x -> (x^k) )

Example 1

This linear second order ODE (Kamke 2.335)is solved by transforming it to the equation for Hypergeometric2F1.

hypergeometricode = 
        4*x*(x  -  1)^2*Derivative[2][y][x]  +  
            2*(x  -  1)*(3*x  -  1)*Derivative[1][y][
                x]  +  (a*x  +  b)*y[x]  ==  0;

sol = Timing[Simplify[y[x] /. DSolve[hypergeometricode, y, x][[1]]]]

Example  2

Here is a third order ODE  (Kamke 3.052) which is solved by a combination of the algorithm for finding exponential solutions, reduction of order and the new equivalence approach to hypergeometric ODEs.

Timing[DSolve[
            x^2 y'''[x]  -  (x  +  3)*x y''[x]  + 
                3*(2*x  +  1) y'[x]  -  
                3*(x  +  1)*y[x]  ==  0,  y,  x]]

Transformations for ODEs with non-rational coefficients

The ODEs which occur in applications (e.g. separation of variables in Laplace's equation) often have trig/hyperbolic/exponential functions in their coefficients.

On the other hand, most modern algorithms (including the ones for Abel ODEs and hypergeometric type ODEs) rely on the coefficients being rational functions.

Mathematica  5.1 includes several techniques for transforming ODEs with non-rational coefficients into ODEs with rational coefficients.

This ODE has a trig function (Cot[x]) in its coefficients and is solved by transforming it to Legendre's equation.

eqn  = y''[x]  +  Cot[x]*y'[x] + 
            v (v + 1)*y[x] == 0;

DSolve[eqn , y[x], x]

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] == 
          0;

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,  Compiled  ->  False,  
      ImageSize  ->  400]

This Riccati equation is solved by transforming it to a linear second order ODE, followed by a change of variable (for Tanh[x]) and finally by appealing to the new code for solving second order hypergeometric ODEs.

Timing[
      sol  =  DSolve[{y'[x]  == y[x]^2  + Tanh[x],  
            y[0] == 1}, y, x]]

Here is a numerical verification of the result.

y[0. ] /. sol[[1]] // Chop

y'[x]  - (y[x]^2  + Tanh[x]) /. 
          sol[[1]] /. {x ->  Random[Complex]} // Chop

Differential equations with Piecewise coefficients

In Mathematica 5.1, DSolve can handle ODEs with UnitStep, Max, Min, etc. (or Piecewise itself) in the coefficients. The equations may have no boundary conditions (general solution) or boundary conditions can be specified to find a particular solution.

Here is the general solution for a linear first order ODE with Piecewise coefficients.

sol = 
      DSolve[ y'[t]  +  
            Piecewise[{{(t^2), t < 0}},  3  t] y[t] == 0,  y,  
        t ]

Plot[y[t] /. sol /. {C[1] -> 4}, {t, -2, 1.5}, 
        PlotStyle -> {Thickness[0.0075], Red}];

Here is the solution for an initial value problem involving different piecewise functions.

sol = 
      DSolve[ y'[x]  +  UnitStep[x - 1]^2 y[x] == 
            Max[x, x^2]  +  Clip[x] &&  y[0] == 1,  y,  
        x ]

Plot[y[x] /. sol, {x, -2, 1.5}, 
        PlotStyle -> {Thickness[0.008], Blue}];

Complete integrals for first order nonlinear PDEs

The general solutions for first order linear and quasilinear PDEs depend on an arbitrary function C[1].

quasilinearPDE  =  
      D[u[x, y], x]  +   D[u[x, y], y] + u[x, y]^2  == 1

Off[Solve::ifun]

DSolve[quasilinearPDE,  u,  {x, y},  
      GeneratedParameters  ->  f]

For nonlinear first order PDEs, there is the notion of a complete integral. This is a solution depending on arbitrary parameters rather than functions.

Mathematica 5.1 includes a variety of techniques for finding complete integrals of nonlinear first order PDEs. The algorithm tries to identify an input equation as belonging to one of the classical types (separable, Clairaut, etc.), perhaps after a transformation (such as Euler and Legendre) and returns the most appropriate solution.

Here is the complete integral for the eikonal equation in two dimensions.

eikonal  =   D[u[x, y], x]^2  +   D[u[x, y], y]^2  == 
        1

sol  =  DSolve[eikonal,  u,  {x, y}][[2]]

The complete integral for the eikonal equation is a two-parameter family of planes.

The following picture illustrates the formation of the envelope of a family of straight lines in the plane. The envelope is a circle.

Plot[
        Evaluate[
          Table[x*Cot[a] + Csc[a],  {a,  1/100, 2 Pi,  
              Pi/16}]],  {x,  -1.4, 1.4},  
        AspectRatio  ->  0.89,  
        PlotRange  ->  {-1.3, 1.3},  
        PlotStyle  ->  {Blue}];

Starting from the complete integral, one can generate new solutions of the PDE by envelope formation i.e.

1. Choose a one-parameter family of solutions (for instance, let C[2] = 5*C[1]).

oneparameterfamily  =  
      u ==  (u[x, y] /. sol[[1]] /.  
            C[2] ->  5*C[1])

2. Form the envelope of this one-parameter family of surfaces.

envelope  =  
      Eliminate [{oneparameterfamily,  
          D[oneparameterfamily,  C[1]]},  {C[1]}]

Then, this envelope is also a solution of the eikonal equation.

The complete integral can be used to solve the characteristic initial value problem for first order nonlinear PDEs.

Summary of new DSolve features in V 5.1

Solutions for equivalence classes of Abel ODEs with non-constant invariants, which can be reduced to inverse linear equations or inverse Riccati equations.

Support for solving linear second order ODEs which are reducible to hypergeometric equations by using a combination of linear transformations in the dependent variable and  bilinear and power transformations in the independent variable.

Transformation techniques for solving ODEs with non-rational (trigonometric, hyperbolic or exponential) coefficients.

Support for finding general solutions and solving initial value problems involving  ODEs with Piecewise  coefficients.

Techniques for finding complete integrals of nonlinear first order PDEs which either belonging to one of the classical types (Clairaut, separable, etc) or which can be reduced to one of these types using Legendre, Euler and other transformations, in the kernel.

Plan for the near future

Provide more detailed documentation for the different aspects of DSolve functionality.

Implement new algorithms for the solution of higher order ODEs as a follow-up to the developments for first and second order ODEs.

Improve the support for solutions of nonlinear PDEs (see the talk by Unal Goktas at this conference).