Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_inteq_abel_weak_weights (d05by)

Purpose

nag_inteq_abel_weak_weights (d05by) computes the fractional quadrature weights associated with the Backward Differentiation Formulae (BDF) of orders 4$4$, 5$5$ and 6$6$. These weights can then be used in the solution of weakly singular equations of Abel type.

Syntax

[wt, sw, ifail] = d05by(iorder, iq, lenfw)
[wt, sw, ifail] = nag_inteq_abel_weak_weights(iorder, iq, lenfw)

Description

nag_inteq_abel_weak_weights (d05by) computes the weights Wi,j${W}_{i,j}$ and ωi${\omega }_{i}$ for a family of quadrature rules related to a BDF method for approximating the integral:
 t 2p − 2 i 1/(sqrt(π)) ∫ (φ(s))/(sqrt(t − s))ds ≃ sqrt(h) ∑ Wi,jφ(j × h) + sqrt(h) ∑ ωi − jφ(j × h),  0 ≤ t ≤ T, 0 j = 0 j = 2p − 1
$1π∫0tϕ(s) t-s ds≃h∑j=0 2p-2Wi,jϕ(j×h)+h∑j=2p-1iωi-jϕ(j×h), 0≤t≤T,$
(1)
with t = i × h(i0)$t=i×h\left(i\ge 0\right)$, for some given h$h$. In (1), p$p$ is the order of the BDF method used and Wi,j${W}_{i,j}$, ωi${\omega }_{i}$ are the fractional starting and the fractional convolution weights respectively. The algorithm for the generation of ωi${\omega }_{i}$ is based on Newton's iteration. Fast Fourier transform (FFT) techniques are used for computing these weights and subsequently Wi,j${W}_{i,j}$ (see Baker and Derakhshan (1987) and Henrici (1979) for practical details and Lubich (1986) for theoretical details). Some special functions can be represented as the fractional integrals of simpler functions and fractional quadratures can be employed for their computation (see Lubich (1986)). A description of how these weights can be used in the solution of weakly singular equations of Abel type is given in Section [Further Comments].

References

Baker C T H and Derakhshan M S (1987) Computational approximations to some power series Approximation Theory (eds L Collatz, G Meinardus and G Nürnberger) 81 11–20
Henrici P (1979) Fast Fourier methods in computational complex analysis SIAM Rev. 21 481–529
Lubich Ch (1986) Discretized fractional calculus SIAM J. Math. Anal. 17 704–719

Parameters

Compulsory Input Parameters

1:     iorder – int64int32nag_int scalar
p$p$, the order of the BDF method to be used.
Constraint: 4iorder6$4\le {\mathbf{iorder}}\le 6$.
2:     iq – int64int32nag_int scalar
Determines the number of weights to be computed. By setting iq to a value, 2iq + 1${2}^{{\mathbf{iq}}+1}$ fractional convolution weights are computed.
Constraint: iq0${\mathbf{iq}}\ge 0$.
3:     lenfw – int64int32nag_int scalar
The dimension of the array wt as declared in the (sub)program from which nag_inteq_abel_weak_weights (d05by) is called.
Constraint: lenfw2iq + 2${\mathbf{lenfw}}\ge {2}^{{\mathbf{iq}}+2}$.

None.

ldsw work lwk

Output Parameters

1:     wt(lenfw) – double array
The first 2iq + 1${2}^{{\mathbf{iq}}+1}$ elements of wt contains the fractional convolution weights ωi${\omega }_{i}$, for i = 0,1,,2iq + 11$\mathit{i}=0,1,\dots ,{2}^{{\mathbf{iq}}+1}-1$. The remainder of the array is used as workspace.
2:     sw(ldsw,2 × iorder1$2×{\mathbf{iorder}}-1$) – double array
ldsw = 2iq + 1 + 2 × iorder1$\mathit{ldsw}={2}^{{\mathbf{iq}}+1}+2×{\mathbf{iorder}}-1$.
sw(i,j + 1)${\mathbf{sw}}\left(\mathit{i},\mathit{j}+1\right)$ contains the fractional starting weights Wi1,j${W}_{\mathit{i}-1,\mathit{j}}$ , for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$ and j = 0,1,,2 × iorder2$\mathit{j}=0,1,\dots ,2×{\mathbf{iorder}}-2$, where n = (2iq + 1 + 2 × iorder1)$\mathit{n}=\left({2}^{{\mathbf{iq}}+1}+2×{\mathbf{iorder}}-1\right)$.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, iorder < 4${\mathbf{iorder}}<4$ or iorder > 6${\mathbf{iorder}}>6$, or iq < 0${\mathbf{iq}}<0$, or lenfw < 2iq + 2${\mathbf{lenfw}}<{2}^{{\mathbf{iq}}+2}$, or ldsw < 2iq + 1 + 2 × iorder − 1$\mathit{ldsw}<{2}^{{\mathbf{iq}}+1}+2×{\mathbf{iorder}}-1$, or lwk < 2iq + 3$\mathit{lwk}<{2}^{{\mathbf{iq}}+3}$.

Accuracy

Not applicable.

Fractional quadrature weights can be used for solving weakly singular integral equations of Abel type. In this section, we propose the following algorithm which you may find useful in solving a linear weakly singular integral equation of the form
 t y(t) = f(t) + 1/(sqrt(π)) ∫ (K(t,s)y(s))/(sqrt(t − s))ds,  0 ≤ t ≤ T, 0
$y(t)=f(t)+1π∫0tK(t,s)y(s) t-s ds, 0≤t≤T,$
(2)
using nag_inteq_abel_weak_weights (d05by). In (2), K(t,s)$K\left(t,s\right)$ and f(t)$f\left(t\right)$ are given and the solution y(t)$y\left(t\right)$ is sought on a uniform mesh of size h$h$ such that T = n × h$T=\mathit{n}×h$. Discretization of (2) yields
 2p − 2 i yi = f(i × h) + sqrt(h) ∑ Wi,jK(i × h,j × h)yj + sqrt(h) ∑ ωi − jK(i × h,j × h)yj, j = 0 j = 2p − 1
$yi = f(i×h) + h ∑ j=0 2p-2 W i,j K (i×h,j×h) yj + h ∑ j=2p-1 i ωi-j K (i×h,j×h) yj ,$
(3)
where yiy(i × h)${y}_{\mathit{i}}\simeq y\left(\mathit{i}×h\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$. We propose the following algorithm for computing yi${y}_{i}$ from (3) after a call to nag_inteq_abel_weak_weights (d05by):
(a) Set n = 2iq + 1 + 2 × iorder2$\mathit{n}={2}^{{\mathbf{iq}}+1}+2×{\mathbf{iorder}}-2$ and h = T / n$h=T/\mathit{n}$.
(b) Equation (3) requires 2 × iorder2$2×{\mathbf{iorder}}-2$ starting values, yj${y}_{\mathit{j}}$, for j = 1,2,,2 × iorder2$\mathit{j}=1,2,\dots ,2×{\mathbf{iorder}}-2$, with y0 = f(0)${y}_{0}=f\left(0\right)$. These starting values can be computed by solving the system
 2 × iorder − 2 yi = f(i × h) + sqrt(h) ∑ sw(i + 1,j + 1)K(i × h,j × h)yj, i = 1,2, … ,2 × iorder − 2. j = 0
$yi = f(i×h) + h ∑ j=0 2×iorder-2 swi+1j+1 K (i×h,j×h) yj , i=1,2,…,2×iorder-2 .$
(c) Compute the inhomogeneous terms
 2 × iorder − 2 σi = f(i × h) + sqrt(h) ∑ sw(i + 1,j + 1)K(i × h,j × h)yj,  i = 2 × iorder − 1,2 × iorder, … ,n. j = 0
$σi = f(i×h) + h ∑ j=0 2×iorder- 2 swi+1j+1 K (i×h,j×h) yj , i = 2 × iorder-1 , 2×iorder , … , n .$
(d) Start the iteration for i = 2 × iorder1,2 × iorder,,n$i=2×{\mathbf{iorder}}-1,2×{\mathbf{iorder}},\dots ,\mathit{n}$ to compute yi${y}_{i}$ from:
 i − 1 (1 − sqrt( h )wt(1)K(i × h,i × h))yi = σi + sqrt(h) ∑ wt(i − j + 1)K(i × h,j × h)yj. j = 2 × iorder − 1
$( 1 - h wt1 K (i×h,i×h) ) yi = σi + h ∑ j=2×iorder-1 i-1 wti-j+1 K (i×h,j×h) yj .$
Note that for nonlinear weakly singular equations, the solution of a nonlinear algebraic system is required at step (b) and a single nonlinear equation at step (d).

Example

```function nag_inteq_abel_weak_weights_example
iorder = int64(4);
iq = int64(3);
lenfw = int64(32);
[wt, sw, ifail] = nag_inteq_abel_weak_weights(iorder, iq, lenfw)
```
```

wt =

0.6928
0.6651
0.4589
0.3175
0.2622
0.2451
0.2323
0.2164
0.2006
0.1878
0.1780
0.1700
0.1629
0.1566
0.1508
0.1457
0.3849
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

sw =

0         0         0         0         0         0         0
0.0565    2.8928   -6.7497   11.6491  -11.1355    5.5374   -1.1223
0.0371    1.7401   -2.8628    6.5207   -6.4058    3.2249   -0.6583
0.0300    1.3207   -2.4642    6.3612   -5.4478    2.7025   -0.5481
0.0258    1.1217   -2.2620    5.3683   -3.7553    2.2132   -0.4549
0.0230    0.9862   -2.0034    4.5005   -3.2772    2.7262   -0.4320
0.0208    0.9001   -1.8989    4.2847   -3.5881    2.8201    0.2253
0.0190    0.8506   -1.9250    4.4164   -4.0181    2.7932    0.1564
0.0173    0.8177   -1.9697    4.5348   -4.2425    2.7458   -0.0697
0.0160    0.7886   -1.9781    4.5318   -4.2769    2.6997   -0.2127
0.0149    0.7603   -1.9548    4.4545   -4.2332    2.6541   -0.2620
0.0140    0.7338   -1.9198    4.3619   -4.1782    2.6059   -0.2716
0.0132    0.7097   -1.8842    4.2754   -4.1246    2.5544   -0.2767
0.0125    0.6880   -1.8497    4.1933   -4.0662    2.5011   -0.2845
0.0119    0.6681   -1.8153    4.1109   -4.0004    2.4479   -0.2915
0.0114    0.6497   -1.7805    4.0279   -3.9304    2.3962   -0.2951
0.0110    0.6327   -1.7461    3.9463   -3.8598    2.3466   -0.2958
0.0105    0.6168   -1.7126    3.8677   -3.7907    2.2990   -0.2950
0.0102    0.6020   -1.6804    3.7926   -3.7238    2.2536   -0.2935
0.0098    0.5882   -1.6495    3.7209   -3.6589    2.2101   -0.2917
0.0095    0.5752   -1.6199    3.6523   -3.5961    2.1686   -0.2895
0.0093    0.5631   -1.5916    3.5867   -3.5356    2.1291   -0.2871
0.0090    0.5517   -1.5644    3.5240   -3.4774    2.0914   -0.2844

ifail =

0

```
```function d05by_example
iorder = int64(4);
iq = int64(3);
lenfw = int64(32);
[wt, sw, ifail] = d05by(iorder, iq, lenfw)
```
```

wt =

0.6928
0.6651
0.4589
0.3175
0.2622
0.2451
0.2323
0.2164
0.2006
0.1878
0.1780
0.1700
0.1629
0.1566
0.1508
0.1457
0.3849
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

sw =

0         0         0         0         0         0         0
0.0565    2.8928   -6.7497   11.6491  -11.1355    5.5374   -1.1223
0.0371    1.7401   -2.8628    6.5207   -6.4058    3.2249   -0.6583
0.0300    1.3207   -2.4642    6.3612   -5.4478    2.7025   -0.5481
0.0258    1.1217   -2.2620    5.3683   -3.7553    2.2132   -0.4549
0.0230    0.9862   -2.0034    4.5005   -3.2772    2.7262   -0.4320
0.0208    0.9001   -1.8989    4.2847   -3.5881    2.8201    0.2253
0.0190    0.8506   -1.9250    4.4164   -4.0181    2.7932    0.1564
0.0173    0.8177   -1.9697    4.5348   -4.2425    2.7458   -0.0697
0.0160    0.7886   -1.9781    4.5318   -4.2769    2.6997   -0.2127
0.0149    0.7603   -1.9548    4.4545   -4.2332    2.6541   -0.2620
0.0140    0.7338   -1.9198    4.3619   -4.1782    2.6059   -0.2716
0.0132    0.7097   -1.8842    4.2754   -4.1246    2.5544   -0.2767
0.0125    0.6880   -1.8497    4.1933   -4.0662    2.5011   -0.2845
0.0119    0.6681   -1.8153    4.1109   -4.0004    2.4479   -0.2915
0.0114    0.6497   -1.7805    4.0279   -3.9304    2.3962   -0.2951
0.0110    0.6327   -1.7461    3.9463   -3.8598    2.3466   -0.2958
0.0105    0.6168   -1.7126    3.8677   -3.7907    2.2990   -0.2950
0.0102    0.6020   -1.6804    3.7926   -3.7238    2.2536   -0.2935
0.0098    0.5882   -1.6495    3.7209   -3.6589    2.2101   -0.2917
0.0095    0.5752   -1.6199    3.6523   -3.5961    2.1686   -0.2895
0.0093    0.5631   -1.5916    3.5867   -3.5356    2.1291   -0.2871
0.0090    0.5517   -1.5644    3.5240   -3.4774    2.0914   -0.2844

ifail =

0

```