bvpFiniteDifference¶

Solves a (parameterized) system of differential equations with boundary conditions at two points, using a variable order, variable step size finite difference method with deferred corrections.

Synopsis¶

bvpFiniteDifference (fcneq, fcnjac, fcnbc, n, nleft, ncupbc, tleft, tright, linear, nfinal, tfinal, yfinal)

Required Arguments¶

void fcneq (n, t, y[], p, dydt[]) (Input): User supplied function to evaluate derivatives.
int n (Input): Number of differential equations.
float t (Input): Independent variable, t.
float y[] (Input): Array of size n containing the dependent variable values, y(t).
float p (Input): Continuation parameter, p. See optional argument problemEmbedded.
float dydt[] (Output): Array of size n containing the derivatives y′(t).
void fcnjac(n, t, y[], p, dypdy[]) (Input): User supplied function to evaluate the Jacobian.
int n (Input): Number of differential equations.
float t (Input): Independent variable, t.
float y[] (Input): Array of size n containing the dependent variable values, y(t).
float p (Input): Continuation parameter, p. See optional argument problemEmbedded.
float dypdy[[]] (Output): n by n array containing the partial derivatives \(a_{i,j} = \partial f_i ⁄ \partial y_j\) evaluated at \((t, y)\). The values \(a_{i,j}\) are returned in dypdy[(i-1)\*n+(j-1)].
void fcnbc(n, yleft[], yright[], p, h[]) (Input): User supplied function to evaluate the boundary conditions.
int n (Input): Number of differential equations.
float yleft[] (Input): Array of size n containing the values of the dependent variable at the left endpoint.
float yright[] (Input): Array of size n containing the values of the dependent variable at the right endpoint.
float p (Input): Continuation parameter, p See optional argument problemEmbedded.
float h[] (Output): Array of size n containing the boundary condition residuals. The boundary conditions are defined by \(h_i = 0\), for i = 0, …, n-1. The left endpoint conditions must be defined first, then, the conditions involving both endpoints, and finally the right endpoint conditions.
int n (Input): Number of differential equations.
int nleft (Input): Number of initial conditions. The value nleft must be greater than or equal to zero and less than n.
int ncupbc (Input): Number of coupled boundary conditions. The value nleft + ncupbc must be greater than zero and less than or equal to n.
float tleft (Input): The left endpoint.
float tright (Input): The right endpoint.
int linear (Input): Integer flag to indicate if the differential equations and the boundary conditions are linear. Set linear to one if the differential equations and the boundary conditions are linear, otherwise set linear to zero.
int nfinal (Output): Number of final grid points, including the endpoints.
float tfinal (Output): Array of size mxgrid containing the final grid points. Only the first nfinal points are significant. See optional argument maxSubinter for definition of mxgrid.
float yfinal (Output): Array of size mxgrid by n containing the values of Y at the points in tfinal. See optional argument maxSubinter for definition of mxgrid.

Optional Arguments¶

tol, float (Input)

Relative error control parameter. The computations stop when

\[\begin{split}\begin{array}{l} \left|E_{i,j}\right| / \max\left(y_{i,j},1.0\right) < \mathit{tol} \text{ for all } i = 0, n = 1 \text{, and } j = 0, \mathit{ngrid}-1 \\ \text{Here } E_{i,j} \text{ is the estimated error on } y_{i,j} \end{array}\end{split}\]

Default: tol = .001.

initialStepsize, int ninit, float tinit[], float yinit[[]], (Input)

Initial gridpoints. Number of initial grid points, including the endpoints, is given by ninit. tinit is an array of size ninit containing the initial grid points. yinit is an array size ninit by n containing an initial guess for the values of Y at the points in tinit.

Default: ninit = 10, tinit[*] equally spaced in the interval [tleft, tright], and yinit[*][*] = 0.

t_print, int (Input)

Parameter indicating the desired output level.

t_print Action

0 No output printed.

1

Intermediate output is printed.

Default: t_print = 0.

maxSubinter, int (Input)

Maximum number of grid points allowed.

Default: mxgrid = 100.

problemEmbedded, float pistep, void fcnpeq(), void fcnpbc()

If this optional argument is supplied, then the function bvpFiniteDifference assumes that the user has embedded the problem into a one-parameter family of problems:

\[h(\mathit{yleft}, \mathit{yright}, p) = 0\]

\[y' = y'(t,y,p)\]

such that for \(p = 0\) the problem is simple. For \(p = 1\), the original problem is recovered. The function bvpFiniteDifference automatically attempts to increment from \(p = 0\) to \(p = 1\). The value pistep is the beginning increment used in this continuation. The increment will usually be changed by function bvpFiniteDifference, but an arbitrary minimum of 0.01 is imposed.

The argument p is the initial increment size for p. The functions fcnpeq and fcnpbc are user-supplied functions, and are defined:

void fcnpeq(n, t, y[], p, dypdp[]) (Input)

User supplied function to evaluate the derivative of y′ with respect to the parameter p.

int n (Input)

Number of differential equations.

float t (Input)

Independent variable, t.

float y[] (Input)

Array of size n containing the dependent variable values.

float p (Input)

Continuation parameter, p.

float dypdp[] (Output)

Array of size n containing the derivative y′ with respect to the parameter p at (t, y).

void fcnpbc(n, yleft[], yright[], p, h[])(Input)

User supplied function to evaluate the derivative of the boundary conditions with respect to the parameter p.

int n (Input)

Number of differential equations.

float yleft[] (Input)

Array of size n containing the values of the dependent variable at the left endpoint.

float yright[] (Input)

Array of size n containing the values of the dependent variable at the right endpoint.

float p (Input)

Continuation parameter, p.

float h[] (Output)

Array of size n containing the derivative of \(f_i\) with respect to p.

errEst (Output)

An array of size n containing estimated error in y.

Description¶

The function bvpFiniteDifference is based on the subprogram PASVA3 by M. Lentini and V. Pereyra (see Pereyra 1978). The basic discretization is the trapezoidal rule over a nonuniform mesh. This mesh is chosen adaptively, to make the local error approximately the same size everywhere. Higher-order discretizations are obtained by deferred corrections. Global error estimates are produced to control the computation. The resulting nonlinear algebraic system is solved by Newton’s method with step control. The linearized system of equations is solved by a special form of Gauss elimination that preserves the sparseness.

Examples¶

Example 1¶

This example solves the third-order linear equation

\[y''' - 2y'' + y' - y = \sin t\]

subject to the boundary conditions \(y(0) = y(2\pi)\) and \(y'(0) = y'(2\pi) = 1\). (Its solution is \(y = \sin t\).) To use bvpFiniteDifference, the problem is reduced to a system of first-order equations by defining \(y_1= y\), \(y_2 = y'\) and \(y_3 = y''\). The resulting system is

\[\begin{split}\begin{array}{ll} y'_1 = y_2 & y_2(0) - 1 = 0 \\ y'_2 = y_3 & y_1(0) - y_1(2\pi) = 0 \\ y'_3 = 2y_3 - y_2 + y_1 + \sin t & y_2(2\pi) - 1 = 0 \\ \end{array}\end{split}\]

Note that there is one boundary condition at the left endpoint \(t = 0\) and one boundary condition coupling the left and right endpoints. The final boundary condition is at the right endpoint. The total number of boundary conditions must be the same as the number of equations (in this case 3).

from __future__ import print_function
from numpy import *
from pyimsl.math.bvpFiniteDifference import bvpFiniteDifference


def fcneqn(n, t, y, p, dydt):
    dydt[0] = y[1]
    dydt[1] = y[2]
    dydt[2] = 2 * y[2] - y[1] + y[0] + sin(t)


def fcnjac(n, t, y, p, dfdy):
    dfdy[0 * n + 0] = 0    # df1/dy1
    dfdy[1 * n + 0] = 0    # df2/dy1
    dfdy[2 * n + 0] = 1    # df3/dy1
    dfdy[0 * n + 1] = 1    # df1/dy2
    dfdy[1 * n + 1] = 0    # df2/dy2
    dfdy[2 * n + 1] = -1   # df3/dy2
    dfdy[0 * n + 2] = 0    # df1/dy3
    dfdy[1 * n + 2] = 1    # df2/dy3
    dfdy[2 * n + 2] = 2    # df3/dy3


def fcnbc(n, yleft, yright, p, h):
    h[0] = yleft[1] - 1
    h[1] = yleft[0] - yright[0]
    h[2] = yright[1] - 1


n = 3
nleft = 1
ncupbc = 1
tleft = 0
tright = 2.0 * pi
linear = 1
yfinal = []
tfinal = []
errest = []
nfinal = []

bvpFiniteDifference(fcneqn, fcnjac, fcnbc,
                    n, nleft, ncupbc, tleft, tright,
                    linear, nfinal, tfinal, yfinal,
                    errEst=errest)

print("          tfinal          y0             y1             y2")
for i in range(0, nfinal[0]):
    print("%5d%15.6e%15.6e%15.6e%15.6e"
          % (i, tfinal[i], yfinal[i][0], yfinal[i][1], yfinal[i][2]))

print("Error estimates     %15.6e%15.6e%15.6e"
      % (errest[0], errest[1], errest[2]))

Output¶

          tfinal          y0             y1             y2
 0.000000e+00  -1.123288e-04   1.000000e+00   6.245272e-05
 3.490659e-01   3.419106e-01   9.397087e-01  -3.419580e-01
 6.981317e-01   6.426907e-01   7.660918e-01  -6.427229e-01
 1.396263e+00   9.847531e-01   1.737334e-01  -9.847452e-01
 2.094395e+00   8.660529e-01  -4.998746e-01  -8.660057e-01
 2.792527e+00   3.421830e-01  -9.395473e-01  -3.420649e-01
 3.490659e+00  -3.417232e-01  -9.396112e-01   3.418946e-01
 4.188790e+00  -8.656880e-01  -5.000588e-01   8.658733e-01
 4.886922e+00  -9.845794e-01   1.734568e-01   9.847519e-01
 5.585054e+00  -6.427722e-01   7.658257e-01   6.429527e-01
 5.934119e+00  -3.420820e-01   9.395433e-01   3.423986e-01
 6.283185e+00  -1.123288e-04   1.000000e+00   6.742588e-04
Error estimates        2.839365e-04   1.792030e-04   5.585375e-04

Example 2¶

In this example, the following nonlinear problem is solved:

\[y'' - y^3 + (1 + sin^2t) sin t = 0\]

with \(y(0) = y(\pi) = 0\). Its solution is \(y = \sin t\). As in Example 1, this equation is reduced to a system of first-order differential equations by defining \(y_1 = y\) and \(y_2 = y'\). The resulting system is

\[\begin{split}\begin{array}{l} y'_1 = y_2y_1(0) = 0 \\ y'_2 = y_1^3 - \left(1+\sin^2t\right) \sin ty_1(\pi) = 0 \end{array}\end{split}\]

In this problem, there is one boundary condition at the left endpoint and one at the right endpoint; there are no coupled boundary conditions.

from __future__ import print_function
from numpy import *
from pyimsl.math.bvpFiniteDifference import bvpFiniteDifference


def fcneqn(n, t, y, p, dydt):
    sx = sin(t)
    dydt[0] = y[1]
    dydt[1] = y[0] * y[0] * y[0] - (sx * sx + 1) * sx


def fcnjac(n, t, y, p, dfdy):
    dfdy[0 * n + 0] = 0            # df1/dy1
    dfdy[1 * n + 0] = 3 * y[0] * y[0]  # df2/dy1
    dfdy[0 * n + 1] = 1            # df1/dy2
    dfdy[1 * n + 1] = 0            # df2/dy2


def fcnbc(n, yleft, yright, p, h):
    h[0] = yleft[0]
    h[1] = yright[0]


n = 2
mxgrid = 100
ninit = 12
nleft = 1
ncupbc = 0
linear = 0
tleft = 0
tright = pi
tinit = zeros((ninit), dtype='double')
yinit = zeros((ninit, n), dtype='double')
tfinal = []
yfinal = []
errest = []
step = (tright - tleft) / (ninit - 1)
nfinal = []

for i in range(0, ninit):
    tinit[i] = tleft + i * step
    yinit[i][0] = 0.4 * (tinit[i] - tleft) * (tright - tinit[i])
    yinit[i][1] = 0.4 * (tright + tleft - 2 * tinit[i])

hinit = {"tinit": tinit, "yinit": yinit}
bvpFiniteDifference(fcneqn, fcnjac, fcnbc,
                    n, nleft, ncupbc, tleft, tright,
                    linear, nfinal, tfinal, yfinal,
                    hinit=hinit, errEst=errest)
print("             t             y0             y1")
for i in range(0, nfinal[0]):
    print("%5d%15.6e%15.6e%15.6e"
          % (i, tfinal[i], yfinal[i][0], yfinal[i][1]))
print("Error estimates     %15.6e%15.6e"
      % (errest[0], errest[1]))

Output¶

             t             y0             y1
 0.000000e+00   0.000000e+00   9.999277e-01
 2.855993e-01   2.817682e-01   9.594315e-01
 5.711987e-01   5.406458e-01   8.412407e-01
 8.567980e-01   7.557380e-01   6.548903e-01
 1.142397e+00   9.096185e-01   4.154530e-01
 1.427997e+00   9.898143e-01   1.423308e-01
 1.713596e+00   9.898143e-01  -1.423308e-01
 1.999195e+00   9.096185e-01  -4.154530e-01
 2.284795e+00   7.557380e-01  -6.548903e-01
 2.570394e+00   5.406458e-01  -8.412407e-01
 2.855993e+00   2.817682e-01  -9.594315e-01
 3.141593e+00   0.000000e+00  -9.999277e-01
Error estimates        3.906276e-05   7.123602e-05

Example 3¶

In this example, the following nonlinear problem is solved:

\[y'' - y^3 = \frac{40}{9}\left(t -\frac{1}{2}\right)^{2/3} - \left(t-\frac{1}{2}\right)^8\]

with \(y(0) = y(1) = \pi/2\). As in the previous examples, this equation is reduced to a system of first-order differential equations by defining \(y_1 = y\) and \(y_2 = y'\). The resulting system is

\[\begin{split}\begin{array}{l} y'_1 = y_2\phantom{...} y_1(0) = \pi / 2 \\ y'_2 = y_1^3 - \tfrac{40}{9}\left(t-\tfrac{1}{2}\right)^{2/3} + \left(t-\tfrac{1}{2}\right)^8 y_1(1) = \pi / 2 \end{array}\end{split}\]

The problem is embedded in a family of problems by introducing the parameter p and by changing the second differential equation to

\[y'_2 = py_1^3 + \frac{40}{9} \left(t-\frac{1}{2}\right)^{2/3} \left(t-\frac{1}{2}\right)^8\]

At \(p = 0\), the problem is linear; and at \(p = 1\), the original problem is recovered. The derivatives \(\partial y' \partial p\) must now be specified in the subroutine fcnpeq. The derivatives \(\partial f / \partial p\) are zero in fcnpbc.

from __future__ import print_function
from numpy import *
from pyimsl.math.bvpFiniteDifference import bvpFiniteDifference


def fcneqn(n, t, y, p, dydt):
    z = t - 0.5
    dydt[0] = y[1]
    dydt[1] = p * y[0] * y[0] * y[0] + 40. / \
        9. * pow(z * z, 1. / 3.) - pow(z, 8)


def fcnjac(n, t, y, p, dfdy):
    dfdy[0 * n + 0] = 0                      # df0/dy0
    dfdy[0 * n + 1] = 1                      # df0/dy1
    dfdy[1 * n + 0] = 3. * (p) * (y[0] * y[0])  # df1/dy0
    dfdy[1 * n + 1] = 0                      # df1/dy1


def fcnbc(n, yleft, yright, p, h):
    pi2 = pi / 2.0
    h[0] = yleft[0] - pi2
    h[1] = yright[0] - pi2


def fcnpeq(n, t, y, p, dfdp):
    dfdp[0] = 0
    dfdp[1] = y[0] * y[0] * y[0]


def fcnpbc(n, yleft, yright, p, dhdp):
    dhdp[0] = 0
    dhdp[1] = 0


n = 2
mxgrid = 45
nleft = 1
ncupbc = 0
tleft = 0
tright = 1
pistep = 0.1
linear = 0
ninit = 5
tinit = [0., 0.4, 0.5, 0.6, 1.0]
yinit = array([[0.15749, 0.00215],
               [0.0, 0.00215],
               [0.15749, -0.83995],
               [-0.05745, 0.0],
               [0.05745, 0.83995]])
tfinal = []
yfinal = []
errest = []
nfinal = []

hinit = {"tinit": tinit, "yinit": yinit}
problem_embedded = {"pistep": pistep, "fcnpeq": fcnpeq,
                    "fcnpbc": fcnpbc}
bvpFiniteDifference(fcneqn, fcnjac, fcnbc, n, nleft,
                    ncupbc, tleft, tright, linear,
                    nfinal, tfinal, yfinal,
                    maxSubinter=mxgrid,
                    problemEmbedded=problem_embedded,
                    hinit=hinit,
                    errEst=errest)
print("             t             y0             y1")
for i in range(0, nfinal[0]):
    print("%5d%15.6e%15.6e%15.6e"
          % (i, tfinal[i], yfinal[i][0], yfinal[i][1]))
print("Error estimates     %15.6e%15.6e"
      % (errest[0], errest[1]))

Output¶

             t             y0             y1
 0.000000e+00   1.570796e+00  -1.949336e+00
 4.444444e-02   1.490495e+00  -1.669566e+00
 8.888889e-02   1.421951e+00  -1.419465e+00
 1.333333e-01   1.363953e+00  -1.194307e+00
 2.000000e-01   1.294526e+00  -8.958462e-01
 2.666667e-01   1.243627e+00  -6.373191e-01
 3.333333e-01   1.208785e+00  -4.135207e-01
 4.000000e-01   1.187783e+00  -2.219352e-01
 4.250000e-01   1.183038e+00  -1.584200e-01
 4.500000e-01   1.179822e+00  -9.973150e-02
 4.625000e-01   1.178748e+00  -7.233901e-02
 4.750000e-01   1.178007e+00  -4.638254e-02
 4.812500e-01   1.177756e+00  -3.399767e-02
 4.875000e-01   1.177582e+00  -2.205550e-02
 4.937500e-01   1.177480e+00  -1.061179e-02
 5.000000e-01   1.177447e+00  -1.566963e-07
 5.062500e-01   1.177480e+00   1.061148e-02
 5.125000e-01   1.177582e+00   2.205519e-02
 5.187500e-01   1.177756e+00   3.399736e-02
 5.250000e-01   1.178007e+00   4.638223e-02
 5.375000e-01   1.178748e+00   7.233870e-02
 5.500000e-01   1.179822e+00   9.973120e-02
 5.750000e-01   1.183038e+00   1.584197e-01
 6.000000e-01   1.187783e+00   2.219349e-01
 6.666667e-01   1.208786e+00   4.135205e-01
 7.333333e-01   1.243628e+00   6.373189e-01
 8.000000e-01   1.294526e+00   8.958459e-01
 8.666667e-01   1.363953e+00   1.194307e+00
 9.111111e-01   1.421951e+00   1.419465e+00
 9.555556e-01   1.490495e+00   1.669566e+00
 1.000000e+00   1.570796e+00   1.949336e+00
Error estimates        3.452906e-06   5.550646e-05

Fatal Errors¶

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”..