plot::ode -- plot the numerical
solution of an ordinary differential equation
Introductionplot::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>], ... )
Parameterst0, 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. |
Optionsmethod |
- | 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. |
Returnsa graphical object of domain type plot::Group.
Related
Functionsnumeric::odesolve, plot, plot::Group, plot::Scene
DetailsY1 := 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.
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.
G.
Cf. example 2.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 = cc
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 = styleG1, 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
1The 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
2We 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
3We 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
4The 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: