Previous Page Next Page Contents

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

Option: Color = c

Option: Style = style

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




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000