d05be:: Integral Equations (NAG Toolbox)

nag_inteq_abel1_weak (d05be) computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the first kind

f (t) + \frac{1}{\sqrt{π}} \int_{0}^{t} \frac{k (t - s)}{\sqrt{t - s}} g (s, y (s)) d s = 0, 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 and if

f (t) = t^{β} w (t) with β > - \frac{1}{2} ​ and ​ w (t) ​ smooth,

(2)

then the solution

y (t)

is unique and has the form

y (t) = t^{β - 1 / 2} z (t)

, (see Lubich (1987)). It is evident from (1) that

f (0) = 0

. You are required to provide the value of

y (t)

t = 0

. If

y (0)

is unknown, 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

4

5

and

6

(

= 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)

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 - 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 + 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_abel1_weak (d05be) is to be called several times (with the same value of

N

) within the same program with different functions.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

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)

is required. This value may be computed from the limit relation (see Gorenflo and Pfeiffer (1991))

\frac{- 2}{\sqrt{π}} k (0) g (0, y (0)) = \lim_{t \to 0} \frac{f (t)}{\sqrt{t}} .

(3)

If the value of the above limit is known then by solving the nonlinear equation (3) an approximation to

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

p

th-order approximation can be used:

\lim_{t \to 0} \frac{f (t)}{\sqrt{t}} ≃ \frac{f (h^{p})}{h^{p / 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

\lim_{t \to 0} \frac{f (t)}{\sqrt{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

2 p - 2

starting values. nag_roots_sys_func_rcomm (c05qd) is used for this purpose. If a failure with

ifail = 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

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

(5)

is required at each step of computation, where

Ψ_{n}

and

α

are constants. nag_inteq_abel1_weak (d05be) calls nag_roots_contfn_cntin_rcomm (c05ax) to find the root of this equation.

Example

The density of the probability that a Brownian motion crosses a one-sided moving boundary

a (t)

before time

t

, satisfies the integral equation (see Hairer et al. (1988))

- \frac{1}{\sqrt{t}} \exp (\frac{1}{2} - {\{a (t)\}}^{2} / t) + \int_{0}^{t} \frac{\exp (- \frac{1}{2} {\{a (t) - a (s)\}}^{2} / (t - s))}{\sqrt{t - s}} y (s) d s = 0, 0 \leq t \leq 7 .

In the case of a straight line

a (t) = 1 + t

, the exact solution is known to be

y (t) = \frac{1}{\sqrt{2 π t^{3}}} \exp \{- {(1 + t)}^{2} / 2 t\}

function d05be_example


fprintf('d05be example results\n\n');

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;
nmesh = 71;
yn = zeros(nmesh,1);
work = zeros(1059, 1);
[yn, work, ifail] = d05be( ...
				 ck, cf, cg, initwt, tlim, yn, work);

h = tlim/double(nmesh-1);
fprintf('Example 1\n\n The stepsize h = %8.4f\n\n',h);
fprintf('     t        Approximate\n');
fprintf('                Solution\n');
t = [0:h:tlim];
sol = [t' yn];
solp = sol(6:5:nmesh,1:2);
fprintf('%8.4f%15.4f\n',solp');

sol1 = @(t)(exp(-(t+1).^2./(2*t))./(sqrt(2*pi)*t.^(3/2)));
[errmax,j] = max(yn - sol1(t)');
fprintf('\nThe maximum absolute error, %10.2e ',errmax);
fprintf(' occurred at t = %8.4f\nwith solution %8.4f\n\n',t(j), yn(j));

ck2 = @(t) sqrt(pi);
cf2 = @(t) -2*log(sqrt(1+t)+sqrt(t))/sqrt(1+t);
cg2 = @(s, y) y;
initwt = 'Subsequent';
tlim = 5;
yn(1) = 1;
[yn, work, ifail] = d05be( ...
			   ck2, cf2, cg2, initwt, tlim, yn, work);
h = tlim/double(nmesh-1);

fprintf('\n\nExample 2\n\n The stepsize h = %8.4f\n\n',h);
fprintf('     t        Approximate\n');
fprintf('                Solution\n');
t = [0:h:tlim];
sol = [t' yn];
solp = sol(8:7:nmesh,1:2);
fprintf('%8.4f%15.4f\n',solp');

sol2 = @(t)(1./(1+t));
[errmax,j] = max(yn - sol2(t)');
fprintf('\nThe maximum absolute error, %10.2e ',errmax);
fprintf(' occurred at t = %8.4f\nwith solution %8.4f\n\n',t(j), yn(j));

d05be example results

Example 1

 The stepsize h =   0.1000

     t        Approximate
                Solution
  0.5000         0.1191
  1.0000         0.0528
  1.5000         0.0265
  2.0000         0.0146
  2.5000         0.0086
  3.0000         0.0052
  3.5000         0.0033
  4.0000         0.0022
  4.5000         0.0014
  5.0000         0.0010
  5.5000         0.0007
  6.0000         0.0004
  6.5000         0.0003
  7.0000         0.0002

The maximum absolute error,   2.86e-03  occurred at t =   0.1000
with solution   0.0326



Example 2

 The stepsize h =   0.0714

     t        Approximate
                Solution
  0.5000         0.6667
  1.0000         0.5000
  1.5000         0.4000
  2.0000         0.3333
  2.5000         0.2857
  3.0000         0.2500
  3.5000         0.2222
  4.0000         0.2000
  4.5000         0.1818
  5.0000         0.1667

The maximum absolute error,   3.17e-06  occurred at t =   0.0714
with solution   0.9333

On entry,	$iorder < 4$ or $iorder > 6$ ,
or	$tlim \leq 10 \times machine precision$ ,
or	$initwt \neq'I'$ or $'S'$ ,
or	$initwt ='S'$ on the first call to nag_inteq_abel1_weak (d05be),
or	$tolnl \leq 10 \times machine precision$ ,
or	$nmesh \neq 2^{m} + 2 \times iorder - 1, m \geq 1$ ,
or	$lwk < (2 \times iorder + 6) \times nmesh + 8 \times {iorder}^{2} - 16 \times iorder + 1$ .

NAG Toolbox: nag_inteq_abel1_weak (d05be)

▸▿ Contents

Purpose

Syntax

Description