NAG Library Function Document

nag_inteq_abel2_weak (d05bdc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

+− 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

nag_inteq_abel2_weak (d05bdc) computes the solution of a weakly singular nonlinear convolution Volterra–Abel integral equation of the second kind using a fractional Backward Differentiation Formulae (BDF) method.

2 Specification

#include <nag.h>

#include <nagd05.h>

void

nag_inteq_abel2_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_abel2_weak (d05bdc) computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the second kind

y (t) = f (t) + \frac{1}{\sqrt{π}} \int_{0}^{t} \frac{k (t - s)}{\sqrt{t - s}} g (s, y (s)) d s, 0 \leq t \leq T .

(1)

Note the constant

\frac{1}{\sqrt{π}}

in (1). It is assumed that the functions involved in (1) are sufficiently smooth.

The function uses a fractional BDF linear multi-step method to generate a family of quadrature rules (see nag_inteq_abel_weak_weights (d05byc)). The BDF methods available in nag_inteq_abel2_weak (d05bdc) are of orders

4

5

and

6

(

= p

say). For a description of the theoretical and practical background to these methods we refer to Lubich (1985) and to Baker and Derakhshan (1987) and Hairer et al. (1988) respectively.

The algorithm is based on computing the solution

y (t)

in a step-by-step fashion on a mesh of equispaced points. The size of the mesh is given by

T / (N - 1)

N

being the number of points at which the solution is sought. These methods require

2 p - 1

(including

y (0)

) 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

N > 32 + 2 p - 1

, 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

N

. An option is provided which avoids the re-evaluation of the fractional weights when nag_inteq_abel2_weak (d05bdc) is to be called several times (with the same value of

N

) within the same program unit 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

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 (1985) Fractional linear multistep methods for Abel–Volterra integral equations of the second kind Math. Comput. 45 463–469

5 Arguments

1: ck – function, supplied by the userExternal Function

ck must evaluate the kernel

k (t)

of the integral equation (1).

The specification of ck is:

double

ck (double t, Nag_Comm *comm)

1: t – doubleInput

On entry:

t

, the value of the independent variable.

2: comm – 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_abel2_weak (d05bdc) you may allocate memory and initialize these pointers with various quantities for use by ck when called from nag_inteq_abel2_weak (d05bdc) (see Section 3.2.1.1 in the Essential Introduction).

2: cf – function, supplied by the userExternal Function

cf must evaluate the function

f (t)

in (1).

The specification of cf is:

double

cf (double t, Nag_Comm *comm)

1: t – doubleInput

On entry:

t

, the value of the independent variable.

2: comm – 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_abel2_weak (d05bdc) you may allocate memory and initialize these pointers with various quantities for use by cf when called from nag_inteq_abel2_weak (d05bdc) (see Section 3.2.1.1 in the Essential Introduction).

3: cg – function, supplied by the userExternal Function

cg must evaluate the function

g (s, y (s))

in (1).

The specification of cg is:

double

cg (double s, double y, Nag_Comm *comm)

1: s – doubleInput

On entry:

s

, the value of the independent variable.

2: y – doubleInput

On entry: the value of the solution

y

at the point s.

3: comm – 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_abel2_weak (d05bdc) you may allocate memory and initialize these pointers with various quantities for use by cg when called from nag_inteq_abel2_weak (d05bdc) (see Section 3.2.1.1 in the Essential Introduction).

4: wtmode – Nag_WeightModeInput

On entry: if the fractional weights required by the method need to be calculated by the function then set

wtmode = Nag_InitWeights

wtmode = Nag_ReuseWeights

, the function assumes the fractional weights have been computed on a previous call and are stored in rwsav.

Constraint:

wtmode = Nag_InitWeights

Nag_ReuseWeights

Note: when nag_inteq_abel2_weak (d05bdc) is re-entered with the value of

wtmode = Nag_ReuseWeights

, the values of nmesh, iorder and the contents of rwsav MUST NOT be changed.

5: iorder – IntegerInput

On entry:

p

, the order of the BDF method to be used.

Suggested value:

iorder = 4

Constraint:

4 \leq iorder \leq 6

6: tlim – doubleInput

On entry: the final point of the integration interval,

T

Constraint:

tlim > 10 \times machine precision

7: tolnl – 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:

tolnl = \sqrt{ε}

where

ε

is the machine precision.

Constraint:

tolnl > 10 \times machine precision

8: nmesh – IntegerInput

On entry:

N

, the number of equispaced points at which the solution is sought.

Constraint:

nmesh = 2^{m} + 2 \times iorder - 1

, where

m \geq 1

9: yn[nmesh] – doubleOutput

On exit:

yn [i - 1]

contains the approximate value of the true solution

y (t)

at the point

t = (i - 1) \times h

, for

i = 1, 2, \dots, nmesh

, where

h = tlim / (nmesh - 1)

10: rwsav[lrwsav] – doubleCommunication Array

On entry: if

wtmode = Nag_ReuseWeights

, rwsav must contain fractional weights computed by a previous call of nag_inteq_abel2_weak (d05bdc) (see description of wtmode).

On exit: contains fractional weights which may be used by a subsequent call of nag_inteq_abel2_weak (d05bdc).

11: lrwsav – IntegerInput

On entry: the dimension of the array rwsav.

Constraint:

lrwsav \geq (2 \times iorder + 6) \times nmesh + 8 \times {iorder}^{2} - 16 \times iorder + 1

12: comm – Nag_Comm *Communication Structure

The NAG communication argument (see Section 3.2.1.1 in the Essential Introduction).

13: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $⟨value⟩$ 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, $iorder = ⟨value⟩$ .
Constraint: $4 \leq iorder \leq 6$ .
NE_INT_2: On entry, $lrwsav = ⟨value⟩$ .
Constraint: $lrwsav \geq (2 \times iorder + 6) \times nmesh + 8 \times {iorder}^{2} - 16 \times iorder + 1$ ; that is, $⟨value⟩$ .
On entry, $nmesh = ⟨value⟩$ and $iorder = ⟨value⟩$ .
Constraint: $nmesh = 2^{m} + 2 \times iorder - 1$ , for some $m$ .
On entry, $nmesh = ⟨value⟩$ and $iorder = ⟨value⟩$ .
Constraint: $nmesh \geq 2 \times iorder + 1$ .
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.
NE_REAL: On entry, $tlim = ⟨value⟩$ .
Constraints: $tlim > 10 \times machine precision$ .
On entry, $tolnl = ⟨value⟩$ .
Constraint: $tolnl > 10 \times machine precision$ .

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 (2) at each step of computation. This value can affect the accuracy of the solution. However, for most problems, the value of

\sqrt{ε}

, where

ε

is the machine precision, should be sufficient.

8 Parallelism and Performance

nag_inteq_abel2_weak (d05bdc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

nag_inteq_abel2_weak (d05bdc) 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 Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

In solving (1), initially, nag_inteq_abel2_weak (d05bdc) computes the solution of a system of nonlinear equations for obtaining the

2 p - 1

starting values. nag_zero_nonlin_eqns_rcomm (c05qdc) is used for this purpose. When a failure with

fail . code =

NE_FAILED_START occurs (which corresponds 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

Y_{n} - α g (t_{n}, Y_{n}) - Ψ_{n} = 0,

(2)

is required at each step of computation, where

Ψ_{n}

and

α

are constants. nag_inteq_abel2_weak (d05bdc) calls nag_zero_cont_func_cntin_rcomm (c05axc) to find the root of this equation.

If a failure with

fail . code =

NE_FAILED_STEP occurs (which corresponds to an error exit from nag_zero_cont_func_cntin_rcomm (c05axc)), you are advised to relax the value of the tolnl or choose a smaller step size by increasing the value of nmesh.

If a failure with

fail . code =

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.