In addition to the default interface to pdsolve/numeric, which uses the default theta methods (the theta family of methods), additional functionality is available for use of numerical pdsolve for educational purposes.

This additional functionality includes the use of 11 standard methods, specification of numerical boundary conditions, specification of startup schemes for two-stage methods, and extensibility for your custom methods.

The use of pdsolve/numeric with these methods is restricted to a single parabolic/hyperbolic PDE that is first order in time.

As with the default method, pdsolve/numeric/method returns a module that can be used to plot solutions, or obtain procedures that provide their values. For information on the plot, animate, $\mathrm{plot3d}$, value, and settings module members, see pdsolve/numeric.

The method argument specifies the method to use for the problem. The available methods are:

For a description of each method, including related problems, and instructions on creating a custom scheme, see pdsolve/numeric/methods.

The options available to pdsolve/numeric/method for these additional methods are as follows.

The time=name, range=l..r, spacestep=numeric, timestep=numeric, bcpts=integer, and optimize=boolean/symbol options are the same as for the default method.

- The time and range options must be specified for problems that are first order in time and space with a boundary condition that involves only one end point (otherwise the time/space variables and range are detected automatically).

- The value of the spacestep option is the size of the spatial step to use for the discrete problem. It defaults to 1/20th the spatial range.

- The value of the timestep option is the size of the temporal step to use for the discrete problem. It defaults to the value of spacestep.

- The value of the bcpts option is the number of points used in the approximation of derivatives present in the boundary conditions.

- The optimize option controls the level of optimization used for generating the solution module. For more information, see pdsolve/numeric.

The numericalbcs=set, list or expression argument describes the numerical boundary conditions. It is required for the use of some methods with some PDEs.

For example, the Euler method relates the values of the solution for 3 consecutive points in the spatial domain. This means that boundary conditions must specify solution values at 2 points in the spatial domain. If the PDE is first order in space, only one natural boundary condition can be specified, and the scheme cannot be used directly.

Numerical boundary conditions allow specification of additional discrete equations, on the spatial mesh, which allow the method to be used for the PDE problem. These are often derived from an extrapolation, or a discretization of the PDE at or near the boundary that requires the additional condition.

The pdsolve command currently accepts two forms of numerical boundary conditions.

The first form specifies that the numerical boundary condition be obtained from one of the specified methods (likely more compact) applied to the input PDE near the boundary.

The second form allows significantly more control, but requires that numerical boundary condition be provided as a discrete equation that is a relationship amongst the discrete values of the dependent variable at the boundary. This form of numerical boundary condition can depend on the time and space variables of the problem, as well as the time and space step sizes (though the names used for the step sizes must be specified), and must be linear in the discrete dependent variable values.

The numericalbcs argument can be specified using the following forms:

where method describes a method accepted by pdsolve (for example, Euler) or a custom method, point describes the boundary and offset to be used for the scheme generated by the method, $\mathrm{timestep}$,$\mathrm{spacestep}$ specify the names used for the time and space step sizes in the numerical boundary condition, and the expression arguments are discrete forms of the numerical boundary conditions.

For the method specified numerical boundary conditions, the value of point indicates the leftmost or rightmost point to be used in the discrete formula obtained by the method. For example, if point is specified as 'n'-1 (where 'n' is some symbol) and the specified method has a width of 2, the condition is applied at the right boundary, and the solution values at the discrete points 'n'-1 and 'n'-2 are used in the discrete formula. As another example, if point is specified as 1 for the same method, the condition is applied at the left boundary, and the solution values at the discrete points 1, 2 are used in the discrete formula.

As a more detailed example, consider solution of the PDE $u_t+0.1u_x=0$ using the LaxFriedrichs method. Since this problem is first order in space, only one natural boundary condition can be applied (one boundary condition on the left boundary). However, the method has a width of 3, so two conditions are needed, thus requiring the specification of a numerical boundary condition.

The additional argument $\mathrm{numericalbcs}=\left[\mathrm{Box}\,n\right]$ specifies the use of the Box method to obtain the required additional relation for the right boundary, and in fact is equivalent to the more complicated specification $\mathrm{numericalbcs}=\left\{\mathrm{spacestep}=h\,\mathrm{timestep}=k\,\frac{v_{1,n}-v_{0,n}+v_{1,n-1}-v_{0,n-1}}{2k}+\frac{0.05000000000\left(v_{1,n}-v_{1,n-1}+v_{0,n}-v_{0,n-1}\right)}{h}=0\right\}$ that corresponds to using expression specified numerical boundary conditions.

The expression specified numerical boundary conditions must be provided in a very particular form. The discrete dependent variable values must be given as indexed variables with the time index first, followed by the space index (that is, $u_{\mathrm{timeidx},\mathrm{spaceidx}}$). Values of the time index can range from $\mathrm{-1}$ to $1$, where a value of $1$ indicates a value at the new time step (Note: The range is from $0$ to $1$ for single stage method). Values of the space index can range from $1$ (left boundary) to $'n'$ (right boundary) where $'n'$ is any symbol that is consistent over the numerical boundary condition. If $'n'$ is present in the space index, it must be of the form $n-i$ where $i\=0,...$ (for example, $u_{1,\frac{n}{2}}$ and $u_{1,2n}$ are not allowed).

For example, a numerical boundary condition that forces the value of the solution on the right boundary to be the same as the value of the solution at the first interior point can be specified as $\mathrm{numericalbcs}=u_{1,n}-u_{1,n-1}$. A similar condition for the left boundary can be specified as $\mathrm{numericalbcs}=u_{1,1}-u_{1,2}$.

The independent variable values also require some description. The spatial values must be given as the space variable indexed in the same manner as the spatial index described for the dependent variable values (for example, $x_{n-1}$ or $x_1$). The temporal values must be given relative to the 0 time index. For example, for the time for the dependent variable value $u_{1,0}$, use the time value $'\mathrm{timevar}'+'\mathrm{timestep}'$ instead of $'\mathrm{timevar}'$.

As a summary example, consider the use of a backward time backward space discretization of the PDE $u_t+u_x-\sin \left(xt\right)$ as a numerical boundary condition for the right boundary. You can obtain the discretizations of the derivatives as $u_t=\frac{u_{1,n}-u_{0,n}}{k}$, and $u_x=\frac{u_{1,n}-u_{1,n-1}}{h}$. Because the expansion point of the discretization is at $t+k,x_n$, replace the time and space variables in the discrete equation with these values. Hence, the specification of this numerical boundary condition is $\mathrm{numericalbcs}=\{\frac{\left(u\[1,n\]-u\[0,n\]\right)}{k}+u_x=\frac{\left(u\[1,n\]-u\[1,n-1\]\right)}{h}-\sin \left(\left(t+k\right)x_n\right),\mathrm{timestep}=k,\mathrm{spacestep}=h\}$ 

Note: This condition can be specified using the method form as $\mathrm{numericalbcs}=\left[{\mathrm{BackwardTime1Space}}_{\mathrm{backward}}\,n\right]$.

Numerical boundary conditions must be provided if needed. Otherwise, Maple returns an error.

The startup=symbol, startupsteps=integer, and startupnumbc=set, list or expression arguments specify how to compute the additional stage required for two stage methods, such as Leapfrog and DuFortFrankel. Methods such as these require the solution values at the current and prior time step to compute the solution at a future time step. This is an issue in the first step of the time integration of the PDE, as only values at one time step are known (that is, the initial condition of the problem).

The startup argument can take two forms, a method specified form, or a function specified form. For the method specified form, the argument must be specified as $\mathrm{startup}=\mathrm{method}$, where method is any of the one stage methods available for use with numerical pdsolve. For the function specified form, the argument must be specified as $\mathrm{startup}=\mathrm{expression}$, where expression is any fully specified function of the time and space variables of the problem. If no startup argument is provided, the method form is used with the default Theta scheme.

The startupsteps argument (default value is $1$) improves the accuracy of the computed solution if the startup method is of lower accuracy than the two stage method. It specifies one greater than the number of doubling steps used to compute the additional stage required by the main method, and is best described through an example.

Consider a PDE where the Leapfrog method is to be used as the primary method, and the Euler method is to be used as the startup method.

* If $\mathrm{startupsteps}=1$, then for the current step size $k$, the Euler method is used to compute the solution value at $t_0+k$, which is then used as the required additional stage.

* If $\mathrm{startupsteps}=2$, then the Euler method is used to compute the solution value at $t_0+\frac{k}{2}$, and the Leapfrog method is used to compute the solution value at $t_0+k$, which is then used as the required additional stage.

* If $\mathrm{startupsteps}=3$, then the Euler method is used to compute the solution value at $t_0+\frac{k}{4}$, and the Leapfrog method is used to compute the solution value at $t_0+\frac{k}{2}$, then at $t_0+k$.

* In general, for $\mathrm{startupsteps}=n$, the Euler method is used to compute the solution value at $t_0+\frac{k}{2^{n-1}}$, then the Leapfrog scheme is used to compute the solution value $n-1$ times at $t_0+\frac{k}{2^{n-2}},t_0+\frac{k}{2^{n-3}},...,t_0+\frac{k}{2},t_0+k$. The value at $t_0+k$ is used as the required additional stage.

This approach allows lower accuracy methods to be used as startup methods with little or no accuracy loss.

The startupnumbc argument specifies the numerical boundary conditions to use for the startup scheme. By default, pdsolve first attempts to obtain the startup solution using the startup method and the same numerical boundary conditions as the primary method. If this fails, pdsolve attempts to obtain the startup solution using the startup method with no numerical boundary conditions. This option is provided to cover cases where neither of these approaches are sufficient. For example, if the primary method requires two numerical boundary conditions, and the startup method requires only one, this approach fails, and startupnumbc must be used. Note that the specification of startupnumbc with the incorrect number of conditions, or use of the option if it is not needed, results in an error.

$\mathrm{PDE1}≔\left\{\mathrm{diff}\left(u\left(x\,t\right)\,t\right)=-\frac{1}{2}\mathrm{diff}\left(u\left(x\,t\right)\,x\right)\right\}$

$\mathrm{PDE1}≔\left\{\frac{\partial }{\partial t}u\left(x\,t\right)=-\frac{\frac{\partial }{\partial x}u\left(x\,t\right)}{2}\right\}$

$\mathrm{IBC1}≔\left\{u\left(0\,t\right)=-\sin \left(\frac{\mathrm{\pi }t}{2}\right)\,u\left(x\,0\right)=\sin \left(\mathrm{\pi }x\right)\right\}$

$\mathrm{IBC1}≔\left\{u\left(0\,t\right)=-\sin \left(\frac{\mathrm{\pi }t}{2}\right)\,u\left(x\,0\right)=\sin \left(\mathrm{\pi }x\right)\right\}$

$\mathrm{smod1}≔\mathrm{pdsolve}\left(\mathrm{PDE1}\,\mathrm{IBC1}\,\mathrm{numeric}\,\mathrm{time}=t\,\mathrm{range}=0..1\,\mathrm{method}=\mathrm{ForwardTime1Space}\left[\mathrm{backward}\right]\right)$

$\mathrm{smod1}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\,\mathrm{plot3d}\,\mathrm{animate}\,\mathrm{value}\,\mathrm{settings}\;...\mathbf{end\; module}$

$\mathrm{esol1}≔\sin \left(\mathrm{\pi }\left(x-\frac{t}{2}\right)\right)$

$\mathrm{esol1}≔\sin \left(\mathrm{\pi }\left(x-\frac{t}{2}\right)\right)$

$\mathrm{errfunc}≔\mathrm{smod1}:-\mathrm{value}\left(u\left(x\,t\right)-\mathrm{esol1}\,t=1.25\,\mathrm{output}=\mathrm{listprocedure}\right)$

$\mathrm{errfunc}≔\left[x=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,t=1.25\,u\left(x\,t\right)-\sin \left(\mathrm{\pi }\left(x-\frac{t}{2}\right)\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\right]$

$\mathrm{errfunc}≔\mathrm{rhs}\left(\mathrm{errfunc}\left[3\right]\right)$

$\mathrm{errfunc}≔\mathbf{proc}\left(x\right)...\mathbf{end\; proc}$

$\mathrm{maxerr}≔\max \left(\mathrm{seq}\left(\mathrm{errfunc}\left(\frac{i}{100}\right)\,i=0..100\right)\right)$

$\mathrm{maxerr}≔0.0299077658092643$

$\mathrm{smod1}:-\mathrm{settings}\left(\mathrm{timestep}=\frac{1}{40}\,\mathrm{spacestep}=\frac{1}{40}\right)\;$$\mathrm{errfunc2}≔\mathrm{smod1}:-\mathrm{value}\left(u\left(x\,t\right)-\mathrm{esol1}\,t=1.25\,\mathrm{output}=\mathrm{listprocedure}\right)$

$\mathrm{errfunc2}≔\left[x=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,t=1.25\,u\left(x\,t\right)-\sin \left(\mathrm{\pi }\left(x-\frac{t}{2}\right)\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\right]$

$\mathrm{errfunc2}≔\mathrm{rhs}\left(\mathrm{errfunc2}\left[3\right]\right)$

$\mathrm{errfunc2}≔\mathbf{proc}\left(x\right)...\mathbf{end\; proc}$

$\mathrm{maxerr2}≔\max \left(\mathrm{seq}\left(\mathrm{errfunc2}\left(\frac{i}{100}\right)\,i=0..100\right)\right)$

$\mathrm{maxerr2}≔0.0156768632175617$

$\frac{\mathrm{maxerr}}{\mathrm{maxerr2}}$

$1.90776467168258$

$\mathrm{PDE2}≔\mathrm{diff}\left(u\left(x\,t\right)\,t\right)+\frac{1}{10}\exp \left(-\frac{t}{10}\right)\mathrm{diff}\left(u\left(x\,t\right)\,x\right)=\frac{\cos \left(\mathrm{\pi }\left(x-\exp \left(-\frac{1}{10}t\right)\right)\right)\mathrm{\pi }\exp \left(-\frac{1}{10}t\right)}{5}$

$\mathrm{PDE2}≔\frac{\partial }{\partial t}u\left(x\,t\right)+\frac{{\ⅇ}^{-\frac{t}{10}}\left(\frac{\partial }{\partial x}u\left(x\,t\right)\right)}{10}=\frac{\cos \left(\mathrm{\pi }\left(x-{\ⅇ}^{-\frac{t}{10}}\right)\right)\mathrm{\pi }{\ⅇ}^{-\frac{t}{10}}}{5}$

$\mathrm{IBC2}≔\left\{u\left(0\,t\right)=u\left(2\,t\right)\,u\left(x\,0\right)=-\sin \left(\mathrm{\pi }x\right)\right\}$

$\mathrm{IBC2}≔\left\{u\left(0\,t\right)=u\left(2\,t\right)\,u\left(x\,0\right)=-\sin \left(\mathrm{\pi }x\right)\right\}$

$\mathrm{discs}≔\mathrm{diff}\left(u\left(x\,t\right)\,t\right)=\frac{u\left[1,n\right]-u\left[0,n\right]+u\left[1,n-1\right]-u\left[0,n-1\right]}{2k},\mathrm{diff}\left(u\left(x\,t\right)\,x\right)=\frac{u\left[1,n\right]-u\left[1,n-1\right]+u\left[0,n\right]-u\left[0,n-1\right]}{2h}$

$\mathrm{discs}≔\frac{\partial }{\partial t}u\left(x\,t\right)=\frac{u_{1,n}-u_{0,n}+u_{1,n-1}-u_{0,n-1}}{2k},\frac{\partial }{\partial x}u\left(x\,t\right)=\frac{u_{1,n}-u_{1,n-1}+u_{0,n}-u_{0,n-1}}{2h}$

$\mathrm{nbc}≔\left\{\mathrm{spacestep}=h\,\mathrm{timestep}=k\,\mathrm{eval}\left(\mathrm{PDE2}\,\left\{\mathrm{discs}\,t=t+\frac{k}{2}\,x=x\left[n\right]-\frac{h}{2}\right\}\right)\right\}$

$\mathrm{nbc}≔\left\{\mathrm{spacestep}=h\,\mathrm{timestep}=k\,\frac{u_{1,n}-u_{0,n}+u_{1,n-1}-u_{0,n-1}}{2k}+\frac{{\ⅇ}^{-\frac{t}{10}-\frac{k}{20}}\left(u_{1,n}-u_{1,n-1}+u_{0,n}-u_{0,n-1}\right)}{20h}=\frac{\cos \left(\mathrm{\pi }\left(x_n-\frac{h}{2}-{\ⅇ}^{-\frac{t}{10}-\frac{k}{20}}\right)\right)\mathrm{\pi }{\ⅇ}^{-\frac{t}{10}-\frac{k}{20}}}{5}\right\}$

$\mathrm{smod2}≔\mathrm{pdsolve}\left(\mathrm{PDE2}\,\mathrm{IBC2}\,\mathrm{type}=\mathrm{numeric}\,\mathrm{numericalbcs}=\mathrm{nbc}\,\mathrm{method}=\mathrm{CrankNicholson}\right)$

$\mathrm{smod2}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\,\mathrm{plot3d}\,\mathrm{animate}\,\mathrm{value}\,\mathrm{settings}\;...\mathbf{end\; module}$

$\mathrm{esol2}≔\sin \left(\mathrm{\pi }\left(x-\exp \left(-\frac{1}{10}t\right)\right)\right)$

$\mathrm{esol2}≔\sin \left(\mathrm{\pi }\left(x-{\ⅇ}^{-\frac{t}{10}}\right)\right)$

$\mathrm{errfunc}≔\mathrm{rhs}\left(\mathrm{smod2}:-\mathrm{value}\left(u\left(x\,t\right)-\mathrm{esol2}\,t=1\,\mathrm{output}=\mathrm{listprocedure}\right)\left[3\right]\right)$

$\mathrm{errfunc}≔\mathbf{proc}\left(x\right)...\mathbf{end\; proc}$

$\mathrm{maxerr}≔\max \left(\mathrm{seq}\left(\mathrm{errfunc}\left(\frac{i}{50}\right)\,i=0..100\right)\right)$

$\mathrm{maxerr}≔0.00486456261189167$

$\mathrm{smod2}:-\mathrm{settings}\left(\mathrm{timestep}=\frac{1}{40}\,\mathrm{spacestep}=\frac{1}{40}\right)\;$$\mathrm{errfunc2}≔\mathrm{rhs}\left(\mathrm{smod2}:-\mathrm{value}\left(u\left(x\,t\right)-\mathrm{esol2}\,t=1\,\mathrm{output}=\mathrm{listprocedure}\right)\left[3\right]\right)$

$\mathrm{errfunc2}≔\mathbf{proc}\left(x\right)...\mathbf{end\; proc}$

$\mathrm{maxerr2}≔\max \left(\mathrm{seq}\left(\mathrm{errfunc2}\left(\frac{i}{50}\right)\,i=0..100\right)\right)$

$\mathrm{maxerr2}≔0.000305388741516921$

$\frac{\mathrm{maxerr}}{\mathrm{maxerr2}}$

$15.9290830032846$

$\mathrm{smod2b}≔\mathrm{pdsolve}\left(\mathrm{PDE2}\,\mathrm{IBC2}\,\mathrm{numeric}\,\mathrm{numericalbcs}=\left[\mathrm{Box}\,n\right]\,\mathrm{method}=\mathrm{CrankNicholson}\right)$

$\mathrm{smod2b}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\,\mathrm{plot3d}\,\mathrm{animate}\,\mathrm{value}\,\mathrm{settings}\;...\mathbf{end\; module}$

$\mathrm{errfunc}≔\mathrm{rhs}\left(\mathrm{smod2b}:-\mathrm{value}\left(u\left(x\,t\right)-\mathrm{esol2}\,t=1\,\mathrm{output}=\mathrm{listprocedure}\right)\left[3\right]\right)$

$\mathrm{errfunc}≔\mathbf{proc}\left(x\right)...\mathbf{end\; proc}$

$\mathrm{maxerr}≔\max \left(\mathrm{seq}\left(\mathrm{errfunc}\left(\frac{i}{50}\right)\,i=0..100\right)\right)$

$\mathrm{maxerr}≔0.00486456261189167$

$\mathrm{PDE3}≔\mathrm{diff}\left(u\left(x\,t\right)\,t\right)=\frac{1}{20}\mathrm{diff}\left(u\left(x\,t\right)\,x\,x\right)+\frac{1}{2}\mathrm{diff}\left(u\left(x\,t\right)\,x\right)$

$\mathrm{PDE3}≔\frac{\partial }{\partial t}u\left(x\,t\right)=\frac{\frac{{\partial }^2}{\partial x^2}u\left(x\,t\right)}{20}+\frac{\frac{\partial }{\partial x}u\left(x\,t\right)}{2}$

$\mathrm{IBC3}≔\left\{u\left(x\,0\right)={\cos \left(\mathrm{\pi }x\right)}^2\,\mathrm{D}\left[1\right]\left(u\right)\left(-\frac{1}{2}\,t\right)=0\,\mathrm{D}\left[1\right]\left(u\right)\left(\frac{1}{2}\,t\right)=0\right\}$

$\mathrm{IBC3}≔\left\{u\left(x\,0\right)={\cos \left(\mathrm{\pi }x\right)}^2\,{\mathrm{D}}_1\left(u\right)\left(-\frac{1}{2}\,t\right)=0\,{\mathrm{D}}_1\left(u\right)\left(\frac{1}{2}\,t\right)=0\right\}$

$\mathrm{smod3}≔\mathrm{pdsolve}\left(\mathrm{PDE3}\,\mathrm{IBC3}\,\mathrm{type}=\mathrm{numeric}\,\mathrm{method}=\mathrm{DuFortFrankel}\,\mathrm{startup}=\mathrm{Euler}\,\mathrm{timestep}=\frac{1}{50}\right)$

$\mathrm{smod3}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\,\mathrm{plot3d}\,\mathrm{animate}\,\mathrm{value}\,\mathrm{settings}\;...\mathbf{end\; module}$

$\mathrm{p1}≔\mathrm{smod3}:-\mathrm{plot}\left(t=0\,\mathrm{color}=\mathrm{red}\right)\:$$\mathrm{p2}≔\mathrm{smod3}:-\mathrm{plot}\left(t=\frac{1}{6}\,\mathrm{color}=\mathrm{maroon}\right)\:$$\mathrm{p3}≔\mathrm{smod3}:-\mathrm{plot}\left(t=\frac{1}{3}\,\mathrm{color}=\mathrm{blue}\right)\:$$\mathrm{p4}≔\mathrm{smod3}:-\mathrm{plot}\left(t=\frac{2}{3}\,\mathrm{color}=\mathrm{green}\right)\:$$\mathrm{p5}≔\mathrm{smod3}:-\mathrm{plot}\left(t=1\,\mathrm{color}=\mathrm{yellow}\right)\:$$\mathrm{plots}\left[\mathrm{display}\right]\left(\left\{\mathrm{p1}\,\mathrm{p2}\,\mathrm{p3}\,\mathrm{p4}\,\mathrm{p5}\right\}\right)$