numeric::odesolve2
--
numerical solution of an ordinary differential equationnumeric::odesolve2
(f, t0, Y0, ..)
returns
a function representing the numerical 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.
numeric::odesolve2(f, t0, Y0 <, method> <, RelativeError = tol> <, Stepsize = h> )
f |
- | procedure representing the vector field of the dynamical system |
t0 |
- | numerical real value for the initial time. |
Y0 |
- | list or 1-dimensional array of numerical values representing the initial value. |
method |
- | name of the numerical scheme, see numeric::odesolve . |
RelativeError =
tol |
- | forces internal numerical Runge-Kutta steps to use
stepsizes with local discretization errors below tol . This
tolerance must be a numerical real value >=
10^(-DIGITS). |
Stepsize = h |
- | switches off the internal error control and uses a
Runge-Kutta iteration with constant stepsize h .
h must be a numerical real value. |
a procedure.
The function returned by numeric::odesolve2
is
sensitive to the environment variable DIGITS
, which determines the
numerical working precision. It uses option remember
.
Y := numeric::odesolve2(f, t0, Y0 <,
options>)
is essentially
Y := t -> numeric::odesolve(t0..t, f, Y0
<,options>).
Y
is called
with a real numerical argument. The call Y(t)
returns the
solution vector in a format corresponding to the type of the initial
condition Y0
with which Y
was defined:
Y(t)
either yields a list or a 1-dimensional array.
If t
is not a real numerical value, then
Y(t)
returns an unevaluated function call.
numeric::odesolve
for details on
the parameters and the options.n
and Symbolic accepted by numeric::odesolve
have no effect:
numeric::odesolve2
ignores these options.The function Y
remembers all values it
has computed. When calling Y(t)
it searches its remember
table for the time T
closest to t
and
integrates from T
to t
using the previously
computed Y(T)
as initial value. This reduces the costs of
a call considerably, if Y
has to be evaluated many times
(e.g., when plotting). However, the result Y(t)
depends on
the history of the MuPAD session! This can lead to unexpected
side effects. Cf. example 3. We recommend to
call Y
only with a monotonically increasing (or
decreasing) sequence of values t
starting from
t0
. Further, the function must be re-initialized whenever
DIGITS
is increased. Cf. example 4.
setuserinfo
(Y,1)
information on the current integration interval is displayed by each
call to Y
.The numerical solution of the initial value problem y'=t*sin(y), y(0)=2 is represented by the following function Y=[y]:
>> f := (t, Y) -> [t*sin(Y[1])]: t0 := 0: Y0 := [2]:
>> Y := numeric::odesolve2(f, 0, [2])
proc Y(t) ... end
It starts numerical integration upon call with a numerical argument:
>> Y(-2), Y(0), Y(0.1), Y(PI + sqrt(2))
[2.968232567], [2.0], [2.004541745], [3.141552691]
Calling Y
with a symbolic argument yields
an unevaluated call:
>> Y(t), Y(t + 5), Y(t^2 - 4)
2 Y(t), Y(t + 5), Y(t - 4)
>> eval(subs(%, t = PI))
[3.13235701], [3.141592654], [3.141592611]
The numerical solution can be plotted. Note that
Y(t)
returns a list, so we plot the list element
Y(t)[1]
:
>> plotfunc2d(Y(t)[1], t = -5..5)
>> delete f, t0, Y0, Y:
We consider the differential equation y''=y^2 with initial conditions y(0)=0, y'(0)=1 The second order equation is converted to a first order system for Y=[Y1,Y2]=[y,y']:
Y1'=Y2, Y2'=Y1^2, Y1(0)=0, Y2(0)=1.
>> f := (t, Y) -> [Y[2], Y[1]^2]: t0 := 0: Y0 := [0, 1]:
>> Y := numeric::odesolve2(f, t0, Y0):
>> Y(1), Y(PI)
[1.087473533, 1.362851121], [1274.867468, 37166.5226]
>> delete f, t0, Y0, Y:
We demonstrate a pitfall with the remember mechanism
built into the functions returned by numeric::odesolve2
.
Consider the initial value problem y'=t*sin(y),
y(0)=2:
>> DIGITS := 5: Y := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [2]):
The following result is stored in the remember table of
Y
:
>> setuserinfo(Y, 1): Y(8.0)
Info: integrating from t0=0 to t=8.0 using Y(t0)=[2] [3.1416]
The following value Y(4.1)
is obtained
using the previously computed Y(8.0)
, integrating backward
in time from t=8.0 to t=4.1:
>> Y(4.1)
Info: integrating from t0=8.0 to t=4.1 using Y(t0)=[3.1416] [4.1634]
We erase the remember table (the fifth operand) of
Y
:
>> Y := subsop(Y, 5 = NIL):
The direct integration from the initial time t0=0 to t=4.1 yields a more accurate result than the previous computation:
>> Y(4.1)
Info: integrating from t0=0 to t=4.1 using Y(t0)=[2] [3.1413]
The reason for this phenomenon becomes obvious, when the solution of the ode is computed for various initial values Y(0)=2,3,4:
>> A := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [2]): B := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [3]): C := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [4]):
>> plotfunc2d(A(t)[1], B(t)[1], C(t)[1], t = 0..8, Grid = 25)
In fact, all solutions with initial values Y(0) in the interval [0,2*PI] approach the same attracting point Y(infinity)=PI. While numerical integration forward in time is a very stable process, it is virtually impossible to recover the correct solution curve integrating backward in time starting at a point close to the attractor.
>> setuserinfo(Y, 0): delete DIGITS, Y, A, B, C:
We consider the system
x'=x+y, y'=x-y, x(0)=1, y(0)=I:
>> f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]: Y := numeric::odesolve2(f, 0, [1, I]):
>> DIGITS := 5: Y(1)
[3.5465 + 1.3683 I, 1.3683 + 0.80989 I]
Increasing DIGITS
does not lead to a more
accurate result because of the remember mechanism:
>> DIGITS := 15: Y(1)
[3.54647799893142 + 1.36829578207979 I, 1.36829578207979 + 0.80988643477182 I]
This is the previous value computed with 5 digits,
printed with 15 digits. Indeed, only 5 digits are correct. For getting
a result that is accurate to full precision one has to erase the
remember table via Y:=subsop(Y,5=NIL)
. Alternatively, one
may create a new numerical solution with a fresh (empty) remember
table:
>> Y := numeric::odesolve2(f, 0, [1, I]): Y(1)
[3.54648242861716 + 1.36829887200859 I, 1.36829887200859 + 0.80988468459998 I]
>> delete f, Y, DIGITS: