hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_inteq_abel1_weak (d05be)

Purpose

nag_inteq_abel1_weak (d05be) 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.

Syntax

[yn, work, ifail] = d05be(ck, cf, cg, initwt, tlim, yn, work, 'iorder', iorder, 'tolnl', tolnl, 'nmesh', nmesh)
[yn, work, ifail] = nag_inteq_abel1_weak(ck, cf, cg, initwt, tlim, yn, work, 'iorder', iorder, 'tolnl', tolnl, 'nmesh', nmesh)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lwk has been removed from the interface
.

Description

nag_inteq_abel1_weak (d05be) computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the first kind
t
f(t) + 1/(sqrt(π))(k(ts))/(sqrt(ts))g(s,y(s))ds = 0,  0tT.
0
f(t)+1π0tk(t-s) t-s g(s,y(s))ds=0,  0tT.
(1)
Note the constant 1/(sqrt(π)) 1π  in (1). It is assumed that the functions involved in (1) are sufficiently smooth and if
f(t) = tβw(t)  with  β > (1/2)​ and ​w(t)​ smooth,
f(t)=tβw(t)  with  β>-12​ and ​w(t)​ smooth,
(2)
then the solution y(t)y(t) is unique and has the form y(t) = tβ1 / 2z(t)y(t)=tβ-1/2z(t), (see Lubich (1987)). It is evident from (1) that f(0) = 0f(0)=0. You are required to provide the value of y(t)y(t) at t = 0t=0. If y(0)y(0) is unknown, Section [Further Comments] 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 (d05by)). The BDF methods available in nag_inteq_abel1_weak (d05be) are of orders 44, 55 and 66 ( = p=p 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 y(t)y(t) in a step-by-step fashion on a mesh of equispaced points. The size of the mesh is given by T / (N1)T/(N-1), NN being the number of points at which the solution is sought. These methods require 2p22p-2 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 + 2p1N>32+2p-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 NN. An option is provided which avoids the re-evaluation of the fractional weights when nag_inteq_abel1_weak (d05be) is to be called several times (with the same value of NN) within the same program with different functions.

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

Parameters

Compulsory Input Parameters

1:     ck – function handle or string containing name of m-file
ck must evaluate the kernel k(t)k(t) of the integral equation (1).
[result] = ck(t)

Input Parameters

1:     t – double scalar
tt, the value of the independent variable.

Output Parameters

1:     result – double scalar
The result of the function.
2:     cf – function handle or string containing name of m-file
cf must evaluate the function f(t)f(t) in (1).
[result] = cf(t)

Input Parameters

1:     t – double scalar
tt, the value of the independent variable.

Output Parameters

1:     result – double scalar
The result of the function.
3:     cg – function handle or string containing name of m-file
cg must evaluate the function g(s,y(s))g(s,y(s)) in (1).
[result] = cg(s, y)

Input Parameters

1:     s – double scalar
ss, the value of the independent variable.
2:     y – double scalar
The value of the solution yy at the point s.

Output Parameters

1:     result – double scalar
The result of the function.
4:     initwt – string (length ≥ 1)
If the fractional weights required by the method need to be calculated by the function then set initwt = 'I'initwt='I' (Initial call).
If initwt = 'S'initwt='S' (Subsequent call), the function assumes the fractional weights have been computed by a previous call and are stored in work.
Constraint: initwt = 'I'initwt='I' or 'S''S'.
Note: when nag_inteq_abel1_weak (d05be) is re-entered with a value of initwt = 'S'initwt='S', the values of nmesh, iorder and the contents of work must not be changed
5:     tlim – double scalar
The final point of the integration interval, TT.
Constraint: tlim > 10 × machine precisiontlim>10×machine precision.
6:     yn(nmesh) – double array
nmesh, the dimension of the array, must satisfy the constraint nmesh = 2m + 2 × iorder1nmesh=2m+2×iorder-1, where m1m1.
yn(1)yn1 must contain the value of y(t)y(t) at t = 0t=0 (see Section [Further Comments]).
7:     work(lwk) – double array
lwk, the dimension of the array, must satisfy the constraint lwk(2 × iorder + 6) × nmesh + 8 × iorder216 × iorder + 1lwk(2×iorder+6)×nmesh+8×iorder2-16×iorder+1.
If initwt = 'S'initwt='S', work must contain fractional weights computed by a previous call of nag_inteq_abel1_weak (d05be) (see description of initwt).

Optional Input Parameters

1:     iorder – int64int32nag_int scalar
pp, the order of the BDF method to be used.
Default: 44
Constraint: 4iorder64iorder6.
2:     tolnl – double scalar
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 [Further Comments]).
Default: sqrt(machine precision)machine precision
Constraint: tolnl > 10 × machine precisiontolnl>10×machine precision.
3:     nmesh – int64int32nag_int scalar
Default: The dimension of the array yn.
NN, the number of equispaced points at which the solution is sought.
Constraint: nmesh = 2m + 2 × iorder1nmesh=2m+2×iorder-1, where m1m1.

Input Parameters Omitted from the MATLAB Interface

lwk nct

Output Parameters

1:     yn(nmesh) – double array
yn(i)yni contains the approximate value of the true solution y(t)y(t) at the point t = (i1) × ht=(i-1)×h, for i = 1,2,,nmeshi=1,2,,nmesh, where h = tlim / (nmesh1)h=tlim/(nmesh-1).
2:     work(lwk) – double array
lwk(2 × iorder + 6) × nmesh + 8 × iorder216 × iorder + 1lwk(2×iorder+6)×nmesh+8×iorder2-16×iorder+1.
Contains fractional weights which may be used by a subsequent call of nag_inteq_abel1_weak (d05be).
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,iorder < 4iorder<4 or iorder > 6iorder>6,
ortlim10 × machine precisiontlim10×machine precision,
orinitwt'I'initwt'I' or 'S''S',
orinitwt = 'S'initwt='S' on the first call to nag_inteq_abel1_weak (d05be),
ortolnl10 × machine precisiontolnl10×machine precision,
ornmesh2m + 2 × iorder1,m1nmesh2m+2×iorder-1,m1,
orlwk < (2 × iorder + 6) × nmesh + 8 × iorder216 × iorder + 1lwk<(2×iorder+6)×nmesh+8×iorder2-16×iorder+1.
  ifail = 2ifail=2
The function cannot compute the 2p22p-2 starting values due to an error in solving the system of nonlinear equations. Relaxing the value of tolnl and/or increasing the value of nmesh may overcome this problem (see Section [Further Comments] for further details).
  ifail = 3ifail=3
The function cannot compute the solution at a specific step due to an error in the solution of the single nonlinear equation (3). Relaxing the value of tolnl and/or increasing the value of nmesh may overcome this problem (see Section [Further Comments] for further details).

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 sqrt(ε)ε, where εε is the machine precision, should be sufficient.

Further Comments

Also when solving (1) the initial value y(0)y(0) is required. This value may be computed from the limit relation (see Gorenflo and Pfeiffer (1991))
(2)/(sqrt(π))k(0)g(0,y(0)) = lim (f(t))/(sqrt(t)).
t0
-2π k(0) g (0,y(0)) = lim t0 f(t) t .
(3)
If the value of the above limit is known then by solving the nonlinear equation (3) an approximation to y(0)y(0) 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 ppth-order approximation can be used:
lim (f(t))/(sqrt(t))(f(hp))/(hp / 2).
t0
lim t0 f(t) t f(hp)hp/2 .
(4)
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 limt0  (f(t))/(sqrt(t)) lim t0 f(t) t  by analytical methods.
Also when solving (1), initially, nag_inteq_abel1_weak (d05be) computes the solution of a system of nonlinear equation for obtaining the 2p22p-2 starting values. nag_roots_sys_func_rcomm (c05qd) is used for this purpose. If a failure with ifail = 2ifail=2 occurs (corresponding to an error exit from nag_roots_sys_func_rcomm (c05qd)), 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
Ynαg(tn,Yn)Ψn = 0,
Yn-αg(tn,Yn)-Ψn=0,
(5)
is required at each step of computation, where ΨnΨn and αα are constants. nag_inteq_abel1_weak (d05be) calls nag_roots_contfn_cntin_rcomm (c05ax) to find the root of this equation.
When a failure with ifail = 3ifail=3 occurs (which corresponds to an error exit from nag_roots_contfn_cntin_rcomm (c05ax)), 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 ifail = 2ifail=2 or 33 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.

Example

function nag_inteq_abel1_weak_example
ck = @(t) exp(-0.5*t);
cf = @(t) -exp( -0.5*(1+t)^2/t)/sqrt(pi*t);
cg = @(s, y) y;
initwt = 'Initial';
tlim = 7;
yn = zeros(71,1);
work = zeros(1059, 1);
[ynOut, workOut, ifail] = nag_inteq_abel1_weak(ck, cf, cg, initwt, tlim, yn, work);
 ynOut, ifail
 

ynOut =

         0
    0.0326
    0.1207
    0.1454
    0.1358
    0.1191
    0.1008
    0.0849
    0.0720
    0.0615
    0.0528
    0.0456
    0.0395
    0.0345
    0.0302
    0.0265
    0.0234
    0.0207
    0.0184
    0.0163
    0.0146
    0.0131
    0.0117
    0.0105
    0.0095
    0.0086
    0.0077
    0.0070
    0.0063
    0.0058
    0.0052
    0.0048
    0.0043
    0.0040
    0.0036
    0.0033
    0.0030
    0.0028
    0.0026
    0.0023
    0.0022
    0.0020
    0.0018
    0.0017
    0.0015
    0.0014
    0.0013
    0.0012
    0.0011
    0.0010
    0.0010
    0.0009
    0.0008
    0.0008
    0.0007
    0.0007
    0.0006
    0.0006
    0.0005
    0.0005
    0.0004
    0.0004
    0.0004
    0.0004
    0.0003
    0.0003
    0.0003
    0.0003
    0.0003
    0.0002
    0.0002


ifail =

                    0


function d05be_example
ck = @(t) exp(-0.5*t);
cf = @(t) -exp( -0.5*(1+t)^2/t)/sqrt(pi*t);
cg = @(s, y) y;
initwt = 'Initial';
tlim = 7;
yn = zeros(71,1);
work = zeros(1059, 1);
[ynOut, workOut, ifail] = d05be(ck, cf, cg, initwt, tlim, yn, work);
 ynOut, ifail
 

ynOut =

         0
    0.0326
    0.1207
    0.1454
    0.1358
    0.1191
    0.1008
    0.0849
    0.0720
    0.0615
    0.0528
    0.0456
    0.0395
    0.0345
    0.0302
    0.0265
    0.0234
    0.0207
    0.0184
    0.0163
    0.0146
    0.0131
    0.0117
    0.0105
    0.0095
    0.0086
    0.0077
    0.0070
    0.0063
    0.0058
    0.0052
    0.0048
    0.0043
    0.0040
    0.0036
    0.0033
    0.0030
    0.0028
    0.0026
    0.0023
    0.0022
    0.0020
    0.0018
    0.0017
    0.0015
    0.0014
    0.0013
    0.0012
    0.0011
    0.0010
    0.0010
    0.0009
    0.0008
    0.0008
    0.0007
    0.0007
    0.0006
    0.0006
    0.0005
    0.0005
    0.0004
    0.0004
    0.0004
    0.0004
    0.0003
    0.0003
    0.0003
    0.0003
    0.0003
    0.0002
    0.0002


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013