NumericalMethods[ODENumeric] - approximates solution of an explicit first order ODE initial problem

Calling Sequence
    ODENumeric(ODE, yinit=y0, x=a..b, method=m, options)
    ODENumeric(ODE, IC, x=a..b, method=m, options)

Parameters
    ODE - equation; a first order differential equation of the form diff(y(x),x)=f(x,y(x))
    y0 - numeric; the value of solution at x=a
    IC - equation; initial condition of the form y(a)=y0
    a..b - range(numeric); the interval where the solution is to be approximated
    m - name; method to be used
    options - (optional) parameters of the form keyword = value, where keyword is one of stepsize, numsteps, output, style, solution, plotops, digits, silent. One of the options stepsize, numsteps must be used.

Description
The ODENumeric command approximates the solution of the initial value problem y' = f(x, y),  y(a) = y0 on the interval [a,b] using classical methods (Euler - forward and implicit, Heun formula, midpoint formula, and the Runge-Kutta method of order 4) and the Taylor method. Its main purpose is to exhibit behaviour of basic methods and be accessible to students not versed in Maple. Acordingly, the main focus is on ease of use and informative output, not performace. Controls are similar to other procedures in the package.

The user must specify the step of the chosen method either directly or by specifying the number of intervals in a partition.

This command is part of the NumericalMethods package, so it can be used in the form ODENumeric(..) only after executing the command with(NumericalMethods). However, it can always be accessed through the long form of the command by using NumericalMethods[ODENumeric](..).

Options

method = foreuler, backeuler, heunform, midpoint, rk4, or [taylor, posint]
- foreuler: the (forward) Euler formula: y_k+1 = y_k + stepsize * f(x_k, y_k).
- backeuler: the implicit (backward) Euler formula:  y_k+1 is determined by solving  y_k+1 = y_k + stepsize * f(x_k+1, y_k+1).
- heunform: the Heun formula (a.k.a. improved Euler formula), a 2nd order Runge-Kutta formula based on trapezoid rule, a predictor-corrector method.
- midpoint: the midpoint method (a.k.a. modified Euler method), a popular Runge-Kutte formula of order 2.
- rk4: the Runge-Kutte formula of order 4.
- [taylor,deg], where deg = posint: the Taylor method. The Taylor polynomial T(t) of degree deg centered at x_k (that is, t = x - x_k ) is constructed using ODE and the initial condition, then y_{k +1} = T(stepsize).

The options argument can be one or more of the following equations.

stepsize = positive
Determines the step size for the iteration. If necessary, the given value is corrected to the nearest step size that can be obtained using (b - a) / n for some integer n. This n then becomes the value for numsteps. If stepsize is not specified, the option numsteps must be used to determine it. If stepsize was either redefined or determined from numsteps, the actual value is reported by the procedure unless prevented by silent.

numsteps = posint
Determines the number of steps (number of intervals in partition) used for the iteration; the step size is then set accordingly and reported by the procedure unless prevented by silent. If numsteps is not specified, the option stepsize must be used.

output = plot, graph, list, information, value, valueerror, or totalerror
Specifies what information is given by the procedure.
- plot returns a plot of the approximate solution.
- graph returns a plot of the approximate solution, and of a quality numeric solution from Maple's numeric solver or the actual solution if supplied by the solution option.
- list returns a list of pairs [x_k , y_k ]. This list can be used directly as argument in the plot command.
- information has an empty return, it prints rows with index k, the value of x_k, the value of the approximation y_k, the value ysol at x_k of Maple's high quality numerical solution or the solution supplied by solution, and the error e_k = | ysol(x_k) - y_k | at x_k .
- value returns the value y_n of the approximate solution at the given endpoint b.
- valueerror returns the error | ysol(b) - y_n |.
- totalerror: returns the global error E_n = max(e_k), where e_k = | ysol(x_k) - y_k |.
The default value is output = plot.

style = line, point, or both
- line plots a graph: a curve consisting of straight segments connecting points [x_k, y_k].
- point plots the points [x_k, y_k].
- both plots both the points [x_k, y_k] and the connecting curve.
The default value is style = line for plots and style = both for graphs.

solution = algebraic
Supplies a function (an expression with variable x) that is to be used for comparison with approximate solution when output is set to graph, information, or totalerror. The procedure checks whether the given expression does solve the given equation. Since deciding equality of two algebraic expressions is not totally reliable, failure of the check does not interrupt the procedure, just a warning is issued unless prevented by silent.

plotops = list
These options are passed to plotting commands if output is set to plot or graph. Some options are overridden by the procedure: the style for approximation is determined by the style option, and if output = graph, then legends are supplied automatically and the solution is shown in line style and color navy.

digits = posint
Specifies how many decimal places should be shown when output = information is specified.
If digits is set to a value greater than Digits, the latter is set to digits for the duration of this procedure. This will influence precision of all calculations during this procedure.
The default value is digits = Digits.

silent = on or off
Controls whether the procedure should print reports while running. Two possible reasons: 1) Solution supplied by solution is incorrect. 2) The stepsize has been set based on numsteps or redefined to fit the range.
- on: such instances are not reported.
- off: such instances are reported.
The default value is silent = off.

Examples

For a quick method introduction, the output = graph option works quite well.
> ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=foreuler, stepsize=0.5, output=graph);

Let's play a bit. Taylor polynomials of degree 1 are tangent lines, so it should again lead to the Euler method and identical plot. We also try the alternative way to supply the initial condition, change step size to force the procedure into adjusting it, change style and try some plot options, the last of which will be ignored.
> ODENumeric(diff(y(x),x)=-y(x), y(0)=3, x=0..5, method=[taylor,1], stepsize=0.513, output=graph, style=line, plotops=[thickness=3,color=gold,style=point]);
Using step size 0.500000.

The point style is useful, as it emphasizes that the actual outcome of a numerical procedure is just points. We will now specify number of steps instead of step size and supply the actual solution.
> ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=backeuler, numsteps=10, output=graph, style=point, solution=3*exp(-x));
Using step size 0.500000.

The default output (plot) is simpler. We will specify an incompatible step size that will be corrected to 0.5, but we supress reporting it with silent.
> ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=backeuler, stepsize=0.505, silent=on);

This simple output allows the user to adjust all its aspects and join it with other graphs easily.
> step:=0.8:x0:=0:B:=4:y0:=3:
  ODEin:=diff(y(x),x)=-y(x), yinit=y0, x=x0..B:
  moreops:=stepsize=step, style=point:
  gsol:=plot(y0*exp(-x), x=x0..B, color=navy, thickness=2,legend="solution"):
  gEf:=ODENumeric(ODEin, method=foreuler, moreops, plotops=[symbol=circle,color=red,symbolsize=18,legend="Euler"]):
  gEb:=ODENumeric(ODEin, method=backeuler, moreops, plotops=[symbol=diamond,color=brown,symbolsize=20,legend="implicit Euler"]):
  gH:=ODENumeric(ODEin, method=heunform, moreops, plotops=[symbol=box,color=blue,symbolsize=15,legend="Heun method (RK2)"]):
  gRK2:=ODENumeric(ODEin, method=midpoint, moreops, plotops=[symbol=cross,color=green,symbolsize=23,legend="midpoint (RK2)"]):
  gRK4:=ODENumeric(ODEin, method=rk4, moreops, plotops=[symbol=circle,color=black,symbolsize=15,legend="RK4"]):
  plots[display]({gsol,gEf,gEb,gH,gRK2,gRK4});

To show the output = information option we will try approximations using cubic Taylor polynomials.
> ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=[taylor,3], stepsize=1, output=information):
k=01:  x= 0.0000000000   y= 3.0000000000   ysol= 3.0000000000   error=0.0000000000
k=02:  x= 1.0000000000   y= 1.0000000000   ysol= 1.1036380554   error=0.1036380554
k=03:  x= 2.0000000000   y= 0.3333333333   ysol= 0.4060056498   error=0.0726723165
k=04:  x= 3.0000000000   y= 0.1111111111   ysol= 0.1493610003   error=0.0382499892
k=05:  x= 4.0000000000   y= 0.0370370370   ysol= 0.0549468538   error=0.0179098168
k=06:  x= 5.0000000000   y= 0.0123456790   ysol= 0.0202138027   error=0.0078681237

It seems that the largest error is 0.1036380554. We will ask for this information directly and supply the actual solution to get precise answer.
> ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=[taylor,3], stepsize=1, output=totalerror, solution=3*exp(-x));

The list output provides data for direct plotting or for further analysis, for instance we can estimate the error of approximation using a control method of higher order.
> data:=ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=foreuler, stepsize=0.50, output=list):
  control:=ODENumeric(diff(y(x),x)=-y(x), yinit=3, x=0..5, method=midpoint, stepsize=0.50, output=list):
  errorest:=[seq([data[i][1],abs(data[i][2]-control[i][2])], i=1..nops(data))]:
  plot(errorest, color=red, thickness=2);

Since we know the analytic solution, it may be interesting to compare this estimate with actual error.
> errorex:=[seq([data[i][1],abs(data[i][2]-3*exp(-data[i][1]))], i=1..nops(data))]:
  plot([errorest,errorex], thickness=2, color=[blue,red], legend=["error estimate","actual error"]);