NAG Library Function Document
nag_inteq_abel1_weak (d05bec)
1
Purpose
nag_inteq_abel1_weak (d05bec) computes the solution of a weakly singular nonlinear convolution Volterra–Abel integral equation of the first kind using a fractional Backward Differentiation Formulae (BDF) method.
2
Specification
#include <nag.h> |
#include <nagd05.h> |
void |
nag_inteq_abel1_weak (
double |
(*ck)(double t,
Nag_Comm *comm),
|
|
double |
(*cf)(double t,
Nag_Comm *comm),
|
|
double |
(*cg)(double s,
double y,
Nag_Comm *comm),
|
|
Nag_WeightMode wtmode,
Integer iorder,
double tlim,
double tolnl,
Integer nmesh,
double yn[],
double rwsav[],
Integer lrwsav,
Nag_Comm *comm,
NagError *fail) |
|
3
Description
nag_inteq_abel1_weak (d05bec) computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the first kind
Note the constant
in
(1). It is assumed that the functions involved in
(1) are sufficiently smooth and if
then the solution
is unique and has the form
, (see
Lubich (1987)). It is evident from
(1) that
. You are required to provide the value of
at
. If
is unknown,
Section 9 gives a description of how an approximate value can be obtained.
The function uses a fractional BDF linear multi-step method selected by you to generate a family of quadrature rules (see
nag_inteq_abel_weak_weights (d05byc)). The BDF methods available in
nag_inteq_abel1_weak (d05bec) are of orders
,
and
(
say). For a description of the theoretical and practical background related to these methods we refer to
Lubich (1987) and to
Baker and Derakhshan (1987) and
Hairer et al. (1988) respectively.
The algorithm is based on computing the solution
in a step-by-step fashion on a mesh of equispaced points. The size of the mesh is given by
,
being the number of points at which the solution is sought. These methods require
starting values which are evaluated internally. The computation of the lag term arising from the discretization of
(1) is performed by fast Fourier transform (FFT) techniques when
, and directly otherwise. The function does not provide an error estimate and you are advised to check the behaviour of the solution with a different value of
. An option is provided which avoids the re-evaluation of the fractional weights when
nag_inteq_abel1_weak (d05bec) is to be called several times (with the same value of
) within the same program with different functions.
4
References
Baker C T H and Derakhshan M S (1987) FFT techniques in the numerical solution of convolution equations J. Comput. Appl. Math. 20 5–24
Gorenflo R and Pfeiffer A (1991) On analysis and discretization of nonlinear Abel integral equations of first kind Acta Math. Vietnam 16 211–262
Hairer E, Lubich Ch and Schlichte M (1988) Fast numerical solution of weakly singular Volterra integral equations J. Comput. Appl. Math. 23 87–98
Lubich Ch (1987) Fractional linear multistep methods for Abel–Volterra integral equations of the first kind IMA J. Numer. Anal 7 97–106
5
Arguments
- 1:
– function, supplied by the userExternal Function
-
ck must evaluate the kernel
of the integral equation
(1).
The specification of
ck is:
double |
ck (double t,
Nag_Comm *comm)
|
|
- 1:
– doubleInput
-
On entry: , the value of the independent variable.
- 2:
– Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
ck.
- user – double *
- iuser – Integer *
- p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_abel1_weak (d05bec) you may allocate memory and initialize these pointers with various quantities for use by
ck when called from
nag_inteq_abel1_weak (d05bec) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: ck should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_abel1_weak (d05bec). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_abel1_weak (d05bec) is likely to produce unexpected results.
- 2:
– function, supplied by the userExternal Function
-
cf must evaluate the function
in
(1).
The specification of
cf is:
double |
cf (double t,
Nag_Comm *comm)
|
|
- 1:
– doubleInput
-
On entry: , the value of the independent variable.
- 2:
– Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cf.
- user – double *
- iuser – Integer *
- p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_abel1_weak (d05bec) you may allocate memory and initialize these pointers with various quantities for use by
cf when called from
nag_inteq_abel1_weak (d05bec) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: cf should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_abel1_weak (d05bec). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_abel1_weak (d05bec) is likely to produce unexpected results.
- 3:
– function, supplied by the userExternal Function
-
cg must evaluate the function
in
(1).
The specification of
cg is:
double |
cg (double s,
double y,
Nag_Comm *comm)
|
|
- 1:
– doubleInput
-
On entry: , the value of the independent variable.
- 2:
– doubleInput
-
On entry: the value of the solution
at the point
s.
- 3:
– Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cg.
- user – double *
- iuser – Integer *
- p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_abel1_weak (d05bec) you may allocate memory and initialize these pointers with various quantities for use by
cg when called from
nag_inteq_abel1_weak (d05bec) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: cg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_abel1_weak (d05bec). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_abel1_weak (d05bec) is likely to produce unexpected results.
- 4:
– Nag_WeightModeInput
-
On entry: if the fractional weights required by the method need to be calculated by the function then set
.
If
, the function assumes the fractional weights have been computed by a previous call and are stored in
rwsav.
Constraint:
or
.
Note: when
nag_inteq_abel1_weak (d05bec) is re-entered with a value of
, the values of
nmesh,
iorder and the contents of
rwsav MUST NOT be changed.
- 5:
– IntegerInput
-
On entry: , the order of the BDF method to be used.
Suggested value:
.
Constraint:
.
- 6:
– doubleInput
-
On entry: the final point of the integration interval, .
Constraint:
.
- 7:
– doubleInput
-
On entry: the accuracy required for the computation of the starting value and the solution of the nonlinear equation at each step of the computation (see
Section 9).
Suggested value:
where is the machine precision.
Constraint:
.
- 8:
– IntegerInput
-
On entry: , the number of equispaced points at which the solution is sought.
Constraint:
, where .
- 9:
– doubleInput/Output
-
On entry:
must contain the value of
at
(see
Section 9).
On exit: contains the approximate value of the true solution at the point , for , where .
- 10:
– doubleCommunication Array
-
On entry: if
,
rwsav must contain fractional weights computed by a previous call of
nag_inteq_abel1_weak (d05bec) (see description of
wtmode).
On exit: contains fractional weights which may be used by a subsequent call of nag_inteq_abel1_weak (d05bec).
- 11:
– IntegerInput
-
On entry: the dimension of the array
rwsav.
Constraint:
.
- 12:
– Nag_Comm *
-
The NAG communication argument (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
- 13:
– NagError *Input/Output
-
The NAG error argument (see
Section 3.7 in How to Use the NAG Library and its Documentation).
6
Error Indicators and Warnings
- NE_ALLOC_FAIL
-
Dynamic memory allocation failed.
See
Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
- NE_BAD_PARAM
-
On entry, argument had an illegal value.
- NE_FAILED_START
-
An error occurred when trying to compute the starting values.
- NE_FAILED_STEP
-
An error occurred when trying to compute the solution at a specific step.
- NE_INT
-
On entry, .
Constraint: .
- NE_INT_2
-
On entry, .
Constraint: ; that is, .
On entry, and .
Constraint: , for some .
On entry, and .
Constraint: .
- NE_INTERNAL_ERROR
-
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
- NE_NO_LICENCE
-
Your licence key may have expired or may not have been installed correctly.
See
Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
- NE_REAL
-
On entry, .
Constraint:
On entry, .
Constraint: .
7
Accuracy
The accuracy depends on
nmesh and
tolnl, the theoretical behaviour of the solution of the integral equation and the interval of integration. The value of
tolnl controls the accuracy required for computing the starting values and the solution of
(3) at each step of computation. This value can affect the accuracy of the solution. However, for most problems, the value of
, where
is the
machine precision, should be sufficient.
8
Parallelism and Performance
nag_inteq_abel1_weak (d05bec) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
nag_inteq_abel1_weak (d05bec) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
x06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
Also when solving
(1) the initial value
is required. This value may be computed from the limit relation (see
Gorenflo and Pfeiffer (1991))
If the value of the above limit is known then by solving the nonlinear equation
(3) an approximation to
can be computed. If the value of the above limit is not known, an approximation should be provided. Following the analysis presented in
Gorenflo and Pfeiffer (1991), the following
th-order approximation can be used:
However, it must be emphasized that the approximation in
(4) may result in an amplification of the rounding errors and hence you are advised (if possible) to determine
by analytical methods.
Also when solving
(1), initially,
nag_inteq_abel1_weak (d05bec) computes the solution of a system of nonlinear equation for obtaining the
starting values.
nag_zero_nonlin_eqns_rcomm (c05qdc) is used for this purpose. If a failure with
NE_FAILED_START occurs (corresponding to an error exit from
nag_zero_nonlin_eqns_rcomm (c05qdc)), you are advised to either relax the value of
tolnl or choose a smaller step size by increasing the value of
nmesh. Once the starting values are computed successfully, the solution of a nonlinear equation of the form
is required at each step of computation, where
and
are constants.
nag_inteq_abel1_weak (d05bec) calls
nag_zero_cont_func_cntin_rcomm (c05axc) to find the root of this equation.
When a failure with
NE_FAILED_STEP occurs (which corresponds to an error exit from
nag_zero_cont_func_cntin_rcomm (c05axc)), you are advised to either relax the value of the
tolnl or choose a smaller step size by increasing the value of
nmesh.
If a failure with
NE_FAILED_START or
NE_FAILED_STEP persists even after adjustments to
tolnl and/or
nmesh then you should consider whether there is a more fundamental difficulty. For example, the problem is ill-posed or the functions in
(1) are not sufficiently smooth.
10
Example
We solve the following integral equations.
Example 1
The density of the probability that a Brownian motion crosses a one-sided moving boundary
before time
, satisfies the integral equation (see
Hairer et al. (1988))
In the case of a straight line
, the exact solution is known to be
Example 2
In this example we consider the equation
The solution is given by
.
In the above examples, the fourth-order BDF is used, and
nmesh is set to
.
10.1
Program Text
Program Text (d05bece.c)
10.2
Program Data
None.
10.3
Program Results
Program Results (d05bece.r)