## `ode` -- the domain of ordinary differential equations

#### Introduction

`ode``(eq, y(x))` represents an ordinary differential equation (ODE) for the function `y(x)`.

`ode``({eq1, eq2, ...}, {y1(x), y2(x), ...})` represents a system of ODEs for the functions `y1(x)`, `y2(x)` etc.

#### Call(s)

`ode(eq, y(x))`
`ode({eq <, inits>}, y(x))`
`ode({eq1, eq2, ... <, inits>}, {y1(x), y2(x), ...})`

#### Parameters

`eq, eq1, eq2, ...` | - | equations or arithmetical expressions in the unknown functions and their derivatives with respect to `x`. An arithmetical expression is regarded as an equation with vanishing right hand side.
`y, y1, y2, ...` | - | the unknown functions: identifiers
`x` | - | the independent variable: an identifier
`inits` | - | the initial or boundary conditions: a sequence of equations

#### Returns

an object of type `ode`.

#### Further Documentation

Section 8.3 of the Tutorial.

#### Related Functions

`numeric::odesolve`, `numeric::odesolve2`, `plot::ode`

#### Details

- In the equations `eq`, `eq1` etc., the unknown functions must be represented by `y(x)`, `y1(x)` etc. Derivatives may be represented either by the `diff` function or by the differential operator `D`. Note that the token `'` provides a handy short cut: `y'(x) = D(y)(x)` = `diff(y(x), x)`.
- The unknown functions must be univariate in the independent variable `x`. Multivariate expressions such as `y(x, t)` are not accepted.
- Initial and boundary conditions are defined by sequences of equations involving the unknown functions or their derivatives on the left hand side. The corresponding values must be specified on the right hand side of the equations. In particular, the differential operator `D` (or the token `'`) must be used to specify values of derivatives at some point. E.g.,

      y(1) = 2, y'(0) = 0, y''(0) = 1

  is a valid sequence of boundary conditions for `inits`.

  Boundary conditions of the first and second kind are allowed. Mixed conditions are not accepted.

  The initial/boundary points and the corresponding initial/boundary values may be symbolic expressions.

- For scalar initial value or boundary value problems, use `ode(`{eq, inits}, y(x)`)` to specify the conditions.
- For systems of ODEs, there must be as many equations as unknown functions.
- The main purpose of the `ode` domain is to provide an environment for overloading the function `solve`.

  In the case of one single equation (possibly together with initial or boundary conditions), `solve` returns a set of explicit solutions or an implicit solution. Each element of the set represents a solution branch.

  In the case of a system of equations, `solve` returns a set of lists of equations for the unknown functions. Each list represents a solution branch.

  An symbolic `solve` call is returned if no solution is found.

- After `setuserinfo(ode, 10)`, a `solve` command provides information on MuPAD's way of solving ODEs.

#### Example 1

In the following, we show how to create and solve a scalar ODE. First, we define the ODE x^2*y'(x) + 3*x*y(x) = sin(x)/x. We use the quote token `'` to represent derivatives:

    >> eq := ode(x^2*y'(x) + 3*x*y(x) = sin(x)/x, y(x))

                    /            sin(x)    2                     \
                 ode| 3 x y(x) - ------ + x  diff(y(x), x), y(x) |
                    \              x                             /

We get an element of the domain `ode` which we can now solve:

    >> solve(eq)

                                 {   C1 + cos(x) }
                                 { - ----------- }
                                 {        3      }
                                 {       x       }

    >> delete eq:

#### Example 2

An initial value problem is defined as a set consisting of the ODE and the initial conditions:

    >> ivp := ode({f''(t) + 4*f(t) = sin(2*t),
                   f(0) = a, f'(0) = b}, f(t))

          ode({f(0) = a, D(f)(0) = b, 4 f(t) + diff(f(t), t, t) -

             sin(2 t)}, f(t))

    >> solve(ivp)

                { sin(2 t)                b sin(2 t)   t cos(2 t) }
                { -------- + a cos(2 t) + ---------- - ---------- }
                {    8                        2            4      }

    >> delete ivp:

#### Example 3

With some restrictions, it is also possible to solve systems of ODEs. First, we define a system:

    >> sys := {x'(t) - x(t) + y(t) = 0, y'(t) - x(t) - y(t) = 0}

              {y(t) - x(t) + D(x)(t) = 0, D(y)(t) - y(t) - x(t) = 0}

A call to `solve` yields the general solution with arbitrary parameters:

    >> solution := solve(ode(sys, {x(t), y(t)}))

          {[y(t) = C7 exp((1 + I) t) + C8 exp((1 - I) t),

             x(t) = I C7 exp((1 + I) t) - I C8 exp((1 - I) t)]}

To verify the result, we substitute it back into the system `sys`. However, for the substitution, it is necessary to rewrite the system into a notation using the `diff` function:

    >> eval(subs(rewrite(sys, diff), op(solution)))

                                      {0 = 0}

    >> delete sys, solution:

#### Example 4

In this example, we point out the various return formats of `ode`'s solve facility. First, we solve an ODE with an initial condition. The solution involves a symbolic integral:

    >> solve(ode({y'(x) + x*y(x) = cos(x), y(0) = 3}, y(x)))

              {      /    2 \                     2 1/2             }
              {      |   x  |   int(cos(t2) exp(t2 )   , t2 = 0..x) }
              { 3 exp| - -- | + ----------------------------------- }
              {      \   2  /                    2 1/2              }
              {                             exp(x )                 }

The following system is solved incompletely:

    >> sys := {x'(t) = -3*y(t)*z(t),
               y'(t) =  3*x(t)*z(t),
               z'(t) = -x(t)*y(t)}:
       solution := solve(ode(sys, {x(t), y(t), z(t)}))

                                      2 1/2                        2 1/2
          {[x(t) = (C10 - C11 + 3 z(t) )   , y(t) = (- C10 - 3 z(t) )   ]

                                          2 1/2
             , [x(t) = (C10 - C11 + 3 z(t) )   ,

                                     2 1/2
             y(t) = - (- C10 - 3 z(t) )   ],

                                          2 1/2
             [x(t) = - (C10 - C11 + 3 z(t) )   ,

                                   2 1/2
             y(t) = (- C10 - 3 z(t) )   ], [

                                         2 1/2
             x(t) = - (C10 - C11 + 3 z(t) )   ,

                                     2 1/2
             y(t) = - (- C10 - 3 z(t) )   ]}

In these four different solutions branches, no solution for the unknown function `z(t)` is provided. In fact, the partial solution above is subject to a further condition on `z(t)`. We substitute the first of our four solutions back into the system `sys`:

    >> eval(subs(rewrite(sys, diff), solution[1]))

          {                                  2 1/2
          { diff(z(t), t) = - (- C10 - 3 z(t) )
          {
          {

                                2 1/2    3 z(t) diff(z(t), t)
             (C10 - C11 + 3 z(t) )   , - -------------------- =
                                                        2 1/2
                                         (- C10 - 3 z(t) )

                                       2 1/2
             3 z(t) (C10 - C11 + 3 z(t) )   ,

               3 z(t) diff(z(t), t)                             2 1/2 }
             ------------------------ = - 3 z(t) (- C10 - 3 z(t) )    }
                                2 1/2                                 }
             (C10 - C11 + 3 z(t) )                                    }

For each solution branch, there remains is exactly one differential equation for `z(t)` to solve.

    >> delete sys, solution:

#### Example 5

It may happen that MuPAD cannot solve a given equation. In such a case, a symbolic `solve` command is returned:

    >> solve(ode(y'(x) + y(x)^2 = b + a*x, y(x)))

                                                          2
                solve(ode(- b - a x + diff(y(x), x) + y(x) , y(x)))

#### Example 6

MuPAD's ODE solver contains algebraic algorithms for computing Liouvillian solutions of linear ordinary differential equations over the rational functions as well. These algorithms are based on differential Galois theory. For second order equations, an algorithm similar to the Kovacic algorithm is implemented. However, instead of computing the so-called Riccati-polynomial of a solution, this algorithm computes the solutions directly using formulas. Only when both solutions are algebraic functions, then in three cases it is necessary to compute the minimal polynomial of a solution, again using formulas. Hence, for Kovacic's famous example

    y''(x) + (3/(16*x^2) + 2/(9*(x - 1)^2) - 3/(16*x*(x - 1)))*y(x) = 0,

it is possible to compute the minimal polynomial of solution:

    >> solve(ode(y''(x) + (3/(16*x^2) + 2/(9*(x - 1)^2)
                           - 3/(16*x*(x - 1)))*y(x), y(x)))

                        24            8        8
          {C2 RootOf(_Y1   + 2985984 x  (x - 1)  -

                      4    8        6         2    16        3
             2799360 x  _Y1  (x - 1)  - 4320 x  _Y1   (x - 1)  -

                     5    4               1/2      1/2 7
             165888 x  _Y1  (x - 2) (I x 3    - I 3   )  +

                     3    12  1/2               1/2      1/2 4
             5760 I x  _Y1   3    (x - 2) (I x 3    - I 3   ) , _Y1)}

MuPAD may find Liouvillian solutions for higher order equations as well. However, there is no guarantee that all of them are found. The following third order equation can be solved completely:

    >> solve(ode(diff(y(x), x, x, x)
               + diff(y(x), x, x)*(2*x^2 - 1)/(2*x^2 - 2*x)
               + diff(y(x), x)*(295*x - 491*x^2 + 196*x^3 - 98)
                                /(196*x^2 - 392*x^3 + 196*x^4)
               + y(x)*(3*x - x^2 - 1)
                       /(196*x^2 - 392*x^3 + 196*x^4), y(x)))

          {
          {                        1/2     1/14
          { C3 (2 x + 2 (x (x - 1))    - 1)     +
          {
          {
          {
          {

                            C4
             -------------------------------- +
                                 1/2     1/14
             (2 x + 2 (x (x - 1))    - 1)

             /
             |                        1/2     1/14
             | C2 (2 x + 2 (x (x - 1))    - 1)
             |
             \

             /
             |                        1/2
             | int(exp(-x) (x (x - 1))
             |
             \

                                 1/2     1/14
             (2 x + 2 (x (x - 1))    - 1)    , x) /

                                 1/2     1/7
             (2 x + 2 (x (x - 1))    - 1)    -

                /                         1/2         \ \ \      }
                |      exp(-x) (x (x - 1))            | | |      }
             int| --------------------------------, x | | | / 28 }
                |                     1/2     1/14    | | |      }
                \ (2 x + 2 (x (x - 1))    - 1)        / / /      }
                                                                 }
                                                                 }

#### Background

- The implemented solution methods mainly stem from
  - Daniel Zwillinger. Handbook of differential equations. San Diego: Academic Press (1992).
  The algebraic algorithms for solving linear ODEs are described in
  - Winfried Fakler. Algebraische Algorithmen zur Lösung von linearen Differentialgleichungen. Stuttgart, Leipzig: Teubner, Reihe MuPAD Reports (1999).
  - Winfried Fakler. On second order homogeneous linear differential equations with Liouvillian solutions. Theor. Comp. Science 187, 27-48 (1997).

#### Changes

- no changes

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000