## `plot::ode` -- plot the numerical solution of an ordinary differential equation

#### Introduction

`plot::ode``([t0, t1, ...], f, Y0, [G])` computes a mesh of numerical sample points Y(t0), Y(t1), ... representing the solution Y(t) of the first order differential equation (dynamical system)

           dY/dt = f(t,Y), Y(t0) = Y0, t0 and t real, Y and Y0 complex vectors.

The procedure

    G: (t, Y) -> [x(t,Y), y(t,Y)] or G: (t, Y) -> [x(t,Y), y(t,Y), z(t,Y)]

maps these solution points (t.i, Y(t.i)) in R x C^n to a mesh of 2D plot points [x.i, y.i] or 3D plot points [x.i, y.i, z.i], respectively,

#### Call(s)

`plot::ode([t0, t1, ...], f, Y0, <method,> <RelativeError = tol,> <Stepsize = h,> [G1 <, Style = style1> <, Color = c1>], [G2 <, Style = style2> <, Color = c2>], ... )`

#### Parameters

`t0, t1, ...` | - | the time mesh: real numerical values. If data are displayed with `Style = Splines`, these values must be in ascending order.
`f` | - | the vector field of the ODE: a procedure. See `numeric::odesolve`.
`Y0` | - | the initial condition of the ODE: a list or a 1-dimensional array. See `numeric::odesolve`.
`G1, G2, ...` | - | ``generators of plot data'': procedures mapping a solution point `(t, Y(t))` to a list `[x, y]` or `[x, y, z]` representing a plot point in 2D or 3D, respectively.

#### Options

`method` | - | use a specific numerical scheme (see `numeric::odesolve`)
`RelativeError = tol` | - | sets a numerical discretization tolerance (see `numeric::odesolve`)
`Stepsize = h` | - | sets a constant stepsize (see `numeric::odesolve`)
`Style = style` | - | sets the style in which the plot data are displayed. The following styles are available: Points, Lines, Splines. The default style is Lines.
`Color = c` | - | sets the RGB color `c` in which the plot data are displayed. The default color is `RGB::Black`.

#### Returns

a graphical object of domain type `plot::Group`.

#### Related Functions

`numeric::odesolve`, `plot`, `plot::Group`, `plot::Scene`

#### Details

- Internally, a sequence of numerical sample points

  `Y1 := numeric::odesolve(t0..t1, f, Y0 <, Options>)`,

  `Y2 := numeric::odesolve(t1..t2, f, Y1 <, Options>)` etc.

  is computed where `Options` is some combination of `method`, `RelativeError= tol`, and `Stepsize= h`. See `numeric::odesolve` for details on the vector field procedure `f`, the initial condition `Y0`, and the options.

- Each of the ``generators of plot data'' `G1`, `G2` etc. creates a graphical solution curve from the sample points `Y0`, `Y1` etc. Each generator `G`, say, is internally called in the form `G(t0, Y0)`, `G(t1, Y1)`, `...` to produce a sequence of plot points in 2D or 3D.

  The solver `numeric::odesolve` returns the solution points `Y0`, `Y1` etc. as lists or 1-dimensional arrays (the actual type is determined by the initial value `Y0`). Consequently, each generator `G` must accept two arguments `(t, Y)`: `t` is a real parameter, `Y` is a ``vector'' (either a list or a 1-dimensional array).

  Each generator must return a list with 2 or 3 elements representing the (x, y) or (x, y, z) coordinates of the graphical point associated with a solution point `(t, Y)` of the ODE. All generators must produce graphical data of the same dimension, i.e., either 2D data as lists with 2 elements, or 3D data as lists with 3 elements.

  Some examples:

  `G := (t, Y) -> [t, Y[1]]` creates a 2D plot of the first component of the solution vector along the y-axis, plotted against the time variable `t` along the x-axis

  `G := (t, Y) -> [Y[1], Y[2]]` creates a 2D phase plot, plotting the first component of the solution along the x-axis and the second component along the y-axis. The result is a solution curve in phase space (parametrized by the time `t`).

  `G := (t, Y) -> [Y[1], Y[2], Y[3]]` creates a 3D phase plot of the first three components of the solution curve.

- Note that arbitrary values associated with the solution curve may be displayed graphically by an appropriate generator `G`. Cf. example 2.
- Several generators `G1`, `G2` etc. can be specified to generate several curves associated with the same numerical mesh `Y0`, `Y1`, `...`. E.g., one can use a generator twice to plot the data with the options `Style = Splines` and `Style = Points`. This allows to display a smooth solution curve together with the sample points. Cf. examples 1, 2, and 3.
- `plot::ode` returns a graphical object of type `plot::Group` containing all graphical solution curves specified by the generators `G1`, `G2` etc. This object can be combined with other plot objects to a graphical scene via `plot::Scene`. Cf. example 3. A final call to `plot` calls the renderer to display the scene.

#### Option: Color = `c`

- This option sets the color in which the graphical data are displayed. The RGB color `c` must be a list of three numerical real values `[r, g, b]` between `0` and `1` representing the red, green and blue contributions according to the RGB color model. A variety of colors is provided by MuPAD's `RGB` data structure.

#### Option: Style = `style`

- The graphical data produced by each of the generators `G1`, `G2` etc. consists of a sequence of mesh points in 2D or 3D, respectively.

  With `Style = Points`, the graphical data are displayed as a discrete set of points.

  With `Style = Lines`, the graphical data points are displayed as a curve consisting of straight line segments between the sample points. The points themselves are not displayed.

  With `Style = Splines`, the graphical data points are displayed as a smooth spline curve connecting the sample points. The points themselves are not displayed.

#### Example 1

The following procedure `f` together with the initial value `Y0` represent the initial value problem dY/dt = f(t, Y) = t*Y - Y^2, Y(t0) = 1. In MuPAD, the function Y is represented by a list with one element:

    >> f := (t, Y) -> [t*Y[1] - Y[1]^2]: Y0 := [1]:

The solution is to be plotted as a function of the time `t`. We define an appropriate generator for the plot data:

    >> G := (t, Y) -> [t, Y[1]]:

The numerical solution is to consist of sample points over the time mesh t.i = i, i = 0,1,..,10. The generator `G` is used twice with different Style options. This generates the sample points together with a smooth spline curve connecting these points:

    >> p := plot::ode([i $ i = 0..10], f, Y0,
                      [G, Style = Points, Color = RGB::Blue],
                      [G, Style = Splines, Color = RGB::Red])

                                   plot::Group()

The point size of the point objects inside `p` is increased:

    >> p:= plot::modify(p, PointWidth = 50):

The object `p` is converted to a graphical scene with various scene options:

    >> s := plot::Scene(p, Labels = ["t", "y"],
              Ticks = [Steps = [1.0, 1], Steps = [1.0, 1]],
              GridLines = Automatic)

                                   plot::Scene()

Finally, the scene is rendered by a call to `plot`:

    >> plot(s)

    >> delete f, Y0, G, p, s:

#### Example 2

We consider the nonlinear oscillator y'' + y + y^3 = sin(t), y(0) = 0, y'(0) = 0.5. As a dynamical system for Y = [y, y'], we have to solve the following initial value problem dY/dt = f(t,Y), Y(0) = Y0:

    >> f := (t, Y) -> [Y[2], sin(t) - Y[1]^3]: Y0 := [0, 0.5]:

The following generator produces a phase plot in the (x,y) plane, embedded in a 3D plot:

    >> G1 := (t, Y) -> [Y[1], Y[2], 0]:

Further, we use the z coordinate of the 3D plot to display the value of the ``energy'' function E = y^2/2 + y'^2/2 over the phase curve:

    >> G2 := (t, Y) -> [Y[1], Y[2], (Y[1]^2 + Y[2]^2)/2]:

The phase curve in the (x, y) plane is combined with the graph of the energy function:

    >> p := plot::ode([i/5 $ i = 0..100], f, Y0,
                      [G1, Style = Splines, Color = RGB::Red],
                      [G2, Style = Points, Color = RGB::Black],
                      [G2, Style = Lines, Color = RGB::Blue]):

We increase the size of the points used in the representation of the energy:

    >> p:= plot::modify(p, PointWidth = 40):

The renderer is called:

    >> plot(p, Axes = Box, Labels = ["y", "y'", "E"],
               CameraPoint = [10, -15, 5]):

    >> delete f, Y0, G1, G2, p:

#### Example 3

We consider the initial value problem y' = f(t, y) = t*sin(t + y^2), y(0) = 1.

    >> f := (t, y) -> t*sin(t + y^2) - y:

The following vector field is tangent to the solution curves:

    >> p1 := plot::vectorfield([1, f(t, y)], t = 1..3, y = 0..1,
                               Grid = [10, 5], Color = RGB::Black):

The following object represents the plot of the solution as a function of `t`:

    >> p2 := plot::ode(
          [i/3 $ i=0..12], (t,Y) -> [f(t, Y[1])], [0],
          [(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Red],
          [(t, Y) -> [t, Y[1]], Style = Splines, Color = RGB::Blue]):

We increase the point size:

    >> p2:= plot::modify(p2, PointWidth = 40):

Finally, we combine the vector field and the ODE plot to a scene and call the renderer:

    >> plot(p1, p2, Scaling = Constrained, Labels = ["t", "y"],
                GridLines = [Steps = 0.25, Steps = 0.25],
                Ticks = [Steps = 1.0, Steps = 0.5]):

    >> delete f, p1, p2:

#### Example 4

The Lorenz ODE is the system

    d/dt [x, y, z] = [p*(y - x), -x*z + r*x - y, x*y - b*z]

with fixed parameters p, r, b. As a dynamical system for Y = [x, y, z], we have to solve the ODE dY/dt = f(t,Y) with the following vector field:

    >> f := proc(t, Y)
            local x, y, z;
            begin
               [x, y, z] := Y:
               [p*(y - x), -x*z + r*x - y, x*y - b*z]
            end_proc:

We consider the following parameters and the following initial condition `Y0`:

    >> p := 10: r := 28: b := 1: Y0 := [1, 1, 1]:

The following generator produces a 3D phase plot:

    >> Gxyz :=  (t, Y) -> Y:

The following generator projects the solution curve to the (y, z) plane with x = -20:

    >> Gyz := (t, Y) -> [-20, Y[2], Y[3]]:

The following generator projects the solution curve to the (x, z) plane with y = 30:

    >> Gxz := (t, Y) -> [Y[1], 30, Y[3]]:

The following generator projects the solution curve to the (x, y) plane with z = 0:

    >> Gxy := (t, Y) -> [Y[1], Y[2], 0]:

With these generators, we create a 3D plot object consisting of the phase curve and its projections. The following call is somewhat time consuming, because it calls the numerical integrator to produce the graphical data:

    >> obj := plot::ode([i/10 $ i=1..100], f, Y0,
               [Gxyz, Style = Splines, Color = RGB::Red],
               [Gyz, Style = Splines, Color = RGB::LightGrey],
               [Gxz, Style = Splines, Color = RGB::LightGrey],
               [Gxy, Style = Splines, Color = RGB::LightGrey]):

Finally, the plot is rendered. Again, this call is somewhat time consuming, because internally, several thousand calls to spline functions occur:

    >> plot(obj, Axes = Box, Ticks = 4,
               CameraPoint = [400, -800, 1200]):

    >> delete f, p, r, b, Y0, Gxyz, Gyz, Gxz, Gxy, obj:

#### Changes

- new function

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