3D plots of ODE solutions
This functionality does not run in MATLAB.
plot::Ode3d(f
,[t_{0}, t_{1}, …]
,Y_{0}
, <[G_{1}, <Style = style_{1}>, <Color = c_{1}>], [G_{2}, <Style = style_{2}>, <Color = c_{2}>], …
>, <method
>, <RelativeError = rtol
>, <AbsoluteError = atol
>, <Stepsize = h
>, <a = a_{min} .. a_{max}
>,options
) plot::Ode3d(f
,[Automatic, t_{start}, t_{end}, t_{step}]
,Y_{0}
, <[G_{1}, <Style = style_{1}>, <Color = c_{1}>], [G_{2}, <Style = style_{2}>, <Color = c_{2}>], …
>, <method
>, <RelativeError = rtol
>, <AbsoluteError = atol
>, <Stepsize = h
>, <a = a_{min} .. a_{max}
>,options
) plot::Ode3d([t_{0}, t_{1}, …]
,f
,Y_{0}
, <[G_{1}, <Style = style_{1}>, <Color = c_{1}>], [G_{2}, <Style = style_{2}>, <Color = c_{2}>], …
>, <method
>, <RelativeError = rtol
>, <AbsoluteError = atol
>, <Stepsize = h
>, <a = a_{min} .. a_{max}
>,options
) plot::Ode3d([Automatic, t_{start}, t_{end}, t_{step}]
,f
,Y_{0}
, <[G_{1}, <Style = style_{1}>, <Color = c_{1}>], [G_{2}, <Style = style_{2}>, <Color = c_{2}>], …
>, <method
>, <RelativeError = rtol
>, <AbsoluteError = atol
>, <Stepsize = h
>, <a = a_{min} .. a_{max}
>,options
)
plot::Ode3d(f, [t_{0}, t_{1},
…], Y_{0})
renders threedimensional
projections of the solutions of the initial value problem given by f
, t_{0}
and Y_{0}
.
plot::Ode3d(f, [t_{0}, t_{1},
…], Y_{0}, [G])
computes a mesh
of numerical sample points Y(t_{0}), Y(t_{1}),
… representing the solution Y(t) of
the first order differential equation (dynamical system)
.
The procedure
maps these solution points (t_{i}, Y(t_{i})) in ℝ×ℂ^{n} to a mesh of 3D plot points [x_{i}, y_{i}, z_{i}]. These points can be connected by straight lines or interpolating splines.
The calling syntax of plot::Ode2d
and plot::Ode3d
as
well as the functionality of these two procedures is identical. The
only difference is that plot::Ode2d
expects graphical
generators G_{1}
, G_{2}
etc.
that produce graphical 2D points, whereas plot::Ode3d
expects
graphical generators producing 3D points.
Internally, a sequence of numerical sample points
Y_1 := numeric::odesolve(f, t_0..t_1, Y_0, Options)
,
Y_2 := numeric::odesolve(f, t_1..t_2, Y_1, Options)
etc.
is computed, where Options
is some combination
of method
, RelativeError = rtol
, AbsoluteError
= atol
, and Stepsize = h
. See numeric::odesolve
for
details on the vector field procedure f
, the initial
condition Y_{0}
, and the options.
The utility function numeric::ode2vectorfield
may
be used to produce the input parameters f, t_{0},
Y_{0}
from a set of differential expressions
representing the ODE. Cf. Example 1.
Each of the "generators of plot data" G_{1}
, G_{2}
etc.
creates a graphical solution curve from the numerical sample points Y_{0}
, Y_{1}
etc.
Each generator G
, say, is internally called in
the form G(t_{0}, Y_{0}),
G(t_{1}, Y_{1}), …
to
produce a sequence of plot points in 3D.
The solver numeric::odesolve
returns the solution
points Y_{0}
, Y_{1}
etc.
as lists or 1dimensional arrays (the actual type is determined by
the initial value Y_{0}
). Consequently,
each generator G
must accept two arguments (t,
Y)
: t
is a real parameter, Y
is
a "vector" (either a list or a 1dimensional array).
Each generator must return a list with 3 elements representing
the (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., for plot::Ode3d
, 3D data as lists with 3
elements.
Some examples:
G := (t, Y) > [Y_1, Y_2, Y_3]
creates
a 3D phase plot of the first three components of the solution curve.
If no generators are given, plot::Ode3d
by
default plots each group of two components as functions of time with
the same style.
Note that arbitrary values associated with the solution curve
may be displayed graphically by an appropriate generator G
.
See Example 2 and Example 5.
Several generators G_{1}, G_{2}
etc.
can be specified to generate several curves associated with the same
numerical mesh Y_{0}, Y_{1},
…
. See Example 1, Example 2,
and Example 3.
The graphical data produced by each of the generators G_{1},
G_{2}
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.
With Style = [Splines, Points]
and Style
= [Lines, Points]
, the effects of the styles used are combined,
i.e., both the evaluation points and the straight lines or splines,
respectively, are displayed.
The plot attributes accepted by plot::Ode2d,Ode3d
include Submesh
= n
, where n is
some positive integer. This attribute only has an effect on the curves
which are returned for the graphical generators with Style
= Splines
and Style = [Splines, Points]
,
respectively. It serves for smoothening the graphical spline curve
using a sufficiently high number of plot points.
n is the number of plot points between two consecutive numerical points corresponding to the time mesh. The default value is n = 4, i.e., the splines are plotted as 5 straight line segments connecting the numerical sample points.
Attribute  Purpose  Default Value 

AbsoluteError  maximal absolute discretization error  
AffectViewingBox  influence of objects on the ViewingBox of
a scene  TRUE 
Colors  list of colors to use  [RGB::Blue , RGB::Red , RGB::Green , RGB::MuPADGold , RGB::Orange , RGB::Cyan , RGB::Magenta , RGB::LimeGreen , RGB::CadmiumYellowLight , RGB::AlizarinCrimson , RGB::Aqua , RGB::Lavender , RGB::SeaGreen , RGB::AureolineYellow , RGB::Banana , RGB::Beige , RGB::YellowGreen , RGB::Wheat , RGB::IndianRed , RGB::Black ] 
Frames  the number of frames in an animation  50 
Function  function expression or procedure  
InitialConditions  initial conditions of the ODE  
Legend  makes a legend entry  
LegendText  short explanatory text for legend  
LegendEntry  add this object to the legend?  FALSE 
LineWidth  width of lines  0.35 
LineStyle  solid, dashed or dotted lines?  Solid 
LinesVisible  visibility of lines  TRUE 
Name  the name of a plot object (for browser and legend)  
ODEMethod  the numerical scheme used for solving the ODE  DOPRI78 
ParameterEnd  end value of the animation parameter  
ParameterName  name of the animation parameter  
ParameterBegin  initial value of the animation parameter  
ParameterRange  range of the animation parameter  
PointSize  the size of points  1.5 
PointStyle  the presentation style of points  FilledCircles 
PointsVisible  visibility of mesh points  TRUE 
Projectors  project an ODE solution to graphical points  
RelativeError  maximal relative discretization error  
Stepsize  set a constant step size  
Submesh  density of submesh (additional sample points)  4 
TimeEnd  end time of the animation  10.0 
TimeMesh  the numerical time mesh  
TimeBegin  start time of the animation  0.0 
TimeRange  the real time span of an animation  0.0 .. 10.0 
Title  object title  
TitleFont  font of object titles  [" sansserif " , 11 ] 
TitlePosition  position of object titles  
TitleAlignment  horizontal alignment of titles w.r.t. their coordinates  Center 
TitlePositionX  position of object titles, x component  
TitlePositionY  position of object titles, y component  
TitlePositionZ  position of object titles, z component  
USubmesh  density of additional sample points for parameter "u"  4 
Visible  visibility  TRUE 
VisibleAfter  object visible after this time value  
VisibleBefore  object visible until this time value  
VisibleFromTo  object visible during this time range  
VisibleAfterEnd  object visible after its animation time ended?  TRUE 
VisibleBeforeBegin  object visible before its animation time starts?  TRUE 
The following procedure f
together with the
initial value Y0
represent the initial value problem
, Y(0)
= 2. In MuPAD^{®}, the 1dimensional vector Y is
represented by a list with one element. The body of the function f
below
addresses the first (and only) entry of this list as Y_{1}
and
returns the 1dimensional vector t Y  Y^{2} as
a list with one element. Also the initial condition Y_{0}
is
a 1dimensional vector represented by a list:
f := (t, Y) > [t*Y[1]  Y[1]^2]: Y0 := [2]:
Alternatively, the utility function numeric::ode2vectorfield
can
be used to generate the input parameters in a more intuitive way:
[f, t0, Y0] := [numeric::ode2vectorfield( {y'(t) = t*y(t)  y(t)^2, y(0) = 2}, [y(t)])]
The numerical solution is to consist of sample points over the
time mesh t_{i} = i, i =
0, 1, …, 10. We use the default generator
of plot::Ode2d
. This generates the sample points
together with a smooth spline curve connecting these points:
p := plot::Ode2d(f, [$ 0..10], Y0, PointSize = 2*unit::mm, PointStyle = Stars):
Finally, the ode solution is rendered by a call to plot
:
plot(p, TicksDistance = 2.0, GridVisible = TRUE, SubgridVisible = TRUE):
We consider the nonlinear oscillator , . As a dynamical system for , we have to solve the following initial value problem , Y(0) = Y_{0}:
f := (t, Y) > [Y[2],  Y[1]^7]: Y0 := [1, 0]:
The following generator produces a plot of the solution Y(t) against the time parameter t:
G1 := (t, Y) > [t, Y[1]]:
Further, we are interested in the values of the function
.
The generator G2
produces the values
along
the solution and plots these values against t:
G2 := (t, Y) > [t, Y[1]^2/2 + Y[2]^2/2]:
The energy function (the "Hamiltonian")
should
be conserved along the solution curve. We define a corresponding generator G3
to
plot
as
a function of t:
G3 := (t, Y) > [t, Y[1]^8/8 + Y[2]^2/2]:
The solution curve is combined with the graph of the function and the conserved energy :
p := plot::Ode2d(f, [i/2 $ i = 0..40], Y0, [G1, Style = Lines, Color = RGB::Red], [G1, Style = Points, Color = RGB::Black], [G2, Style = Lines, Color = RGB::Blue], [G2, Style = Points, Color = RGB::Black], [G3, Style = Lines, Color = RGB::Green], [G3, Style = Points, Color = RGB::Black], PointSize = 1.5*unit::mm, LineWidth = 0.2*unit::mm ):
Note that by using each generator twice, we are able to set different colors for the lines and points. The renderer is called:
plot(p):
To visualize the dependency of the trajectory on the initial
conditions, we animate plot::Ode2d
over different
values of
:
plot(plot::Ode2d(f, [i/6 $ i = 0..120], [1, a], a = 1/2..1/2, [G1, Style = Lines, Color = RGB::Red], [G2, Style = Lines, Color = RGB::Blue], [G3, Style = Lines, Color = RGB::Green], LineWidth = 0.2*unit::mm, Frames=25))
We consider the initial value problem , y(0) = 0:
f := (t, y) > t*sin(t + y^2): Y0:= [0]:
The following vector field is tangent to the solution curves:
p1 := plot::VectorField2d([1, f(t, y)], t = 0..4, y = 1.2..1.2, Mesh = [21, 25], Color = RGB::Black):
The following object represents the plot of the solution as
a function of t
:
p2 := plot::Ode2d( (t,Y) > [f(t, Y[1])], [i/3 $ i=0..12], Y0, [(t, Y) > [t, Y[1]], Style = Points, Color = RGB::Red], [(t, Y) > [t, Y[1]], Style = Splines, Color = RGB::Blue]):
We define the point size explicitly:
p2::PointSize := 2*unit::mm:
Finally, we combine the vector field and the ODE plot to a scene and call the renderer:
plot(p1, p2, XTicksDistance = 0.5, YTicksDistance = 0.2, Axes = Frame, AxesTitles = ["t", "y"], GridVisible = TRUE):
By default, numeric::odesolve
(which is used by plot::Ode2d
and plot::Ode3d
internally)
uses adaptive step sizes and a method of order 8. Usually, there is
no reason to change these settings, except for demonstrative purposes.
In the following animation, we use a straightforward explicit Euler
method (of first order) and show how decreasing the step size improves
the quality of the calculated solution.
Our differential equation is , obviously fulfilled by the exponential function:
[f, t0, Y0] := [numeric::ode2vectorfield( {y'(t)=y(t), y(0)=1}, [y(t)])]:
To judge the quality of the numerical solution, we plot the symbolic solution alongside the approximation:
plot(plot::Function2d(exp(x), x=0..3, Color = RGB::Black, LineStyle = Dashed), plot::Ode2d(f, [Automatic, 0, 3, 1/n], Y0, n = 1..50, EULER1, Stepsize = 1/n, [(t, Y) > [t, Y[1]], Style=[Lines, Points]]))
We consider the nonlinear oscillator , . As a dynamical system for , we have to solve the following initial value problem , Y(0) = Y_{0}:
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 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::Ode3d(f, [i/5 $ i = 0..100], Y0, [G1, Style = Splines, Color = RGB::Red], [G2, Style = Points, Color = RGB::Black], [G2, Style = Lines, Color = RGB::Blue]):
We set an explicit size of the points used in the representation of the energy:
p::PointSize := 2*unit::mm:
The renderer is called:
plot(p, AxesTitles = ["y", "y'", "E"], CameraDirection = [10, 15, 5]):
The Lorenz ODE is the system
with fixed parameters p, r, b. As a dynamical system for Y = [x, y, z], we have to solve the ODE 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 Gxyz
produces a 3D
phase plot of the solution. The generator Gyz
projects
the solution curve to the (y, z) plane
with x = 20;
the generator Gxz
projects the solution curve to
the (x, z) plane
with y =  15;
the generator Gxy
projects the solution curve to
the (x, y) plane
with z = 0:
Gxyz := (t, Y) > Y: Gyz := (t, Y) > [ 20, Y[2], Y[3]]: Gxz := (t, Y) > [Y[1], 15, Y[3]]: 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.
object := plot::Ode3d(f, [i/10 $ i=1..100], Y0, [Gxyz, Style = Splines, Color = RGB::Red], [Gyz, Style = Splines, Color = RGB::Grey50], [Gxz, Style = Splines, Color = RGB::Grey50], [Gxy, Style = Splines, Color = RGB::Grey50], Submesh = 7):
Finally, the plot is rendered. This call is somewhat time consuming
because it calls the numerical solver numeric::odesolve
to
produce the graphical data:
plot(object, CameraDirection = [220, 110, 150])

The vector field of the ODE: a procedure.
See


The time mesh: real numerical values. If data are displayed
with


The time mesh: real numerical values.


The initial condition of the ODE: a list or
a 1dimensional array. See


"generators of plot data": procedures mapping
a solution point


Use a specific numerical scheme (see 

Animation parameter, specified as 

Option, specified as Sets the style in which the plot data are displayed. The following
styles are available: 

Option, specified as Sets the RGB
color 

Option, specified as Sets a numerical discretization tolerance (see 

Option, specified as Sets a numerical discretization tolerance (see 

Option, specified as Sets a constant stepsize (see 