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_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 44, 55 and 66. 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,jWi,j and ωiωi for a family of quadrature rules related to a BDF method for approximating the integral:
t 2p2 i
1/(sqrt(π))(φ(s))/(sqrt(ts))dssqrt(h)Wi,jφ(j × h) + sqrt(h)ωijφ(j × h),  0tT,
0 j = 0 j = 2p1
1π0tϕ(s) t-s dshj=0 2p-2Wi,jϕ(j×h)+hj=2p-1iωi-jϕ(j×h),  0tT,
(1)
with t = i × h(i0)t=i×h(i0), for some given hh. In (1), pp is the order of the BDF method used and Wi,jWi,j, ωiωi are the fractional starting and the fractional convolution weights respectively. The algorithm for the generation of ωiωi is based on Newton's iteration. Fast Fourier transform (FFT) techniques are used for computing these weights and subsequently Wi,jWi,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
pp, the order of the BDF method to be used.
Constraint: 4iorder64iorder6.
2:     iq – int64int32nag_int scalar
Determines the number of weights to be computed. By setting iq to a value, 2iq + 12iq+1 fractional convolution weights are computed.
Constraint: iq0iq0.
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 + 2lenfw2iq+2.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

ldsw work lwk

Output Parameters

1:     wt(lenfw) – double array
The first 2iq + 12iq+1 elements of wt contains the fractional convolution weights ωiωi, for i = 0,1,,2iq + 11i=0,1,,2iq+1-1. The remainder of the array is used as workspace.
2:     sw(ldsw,2 × iorder12×iorder-1) – double array
ldsw = 2iq + 1 + 2 × iorder1ldsw=2iq+1+2×iorder-1.
sw(i,j + 1)swij+1 contains the fractional starting weights Wi1,jWi-1,j , for i = 1,2,,ni=1,2,,n and j = 0,1,,2 × iorder2j=0,1,,2×iorder-2, where n = (2iq + 1 + 2 × iorder1)n=(2iq+1+2×iorder-1).
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,
oriq < 0iq<0,
orlenfw < 2iq + 2lenfw<2iq+2,
orldsw < 2iq + 1 + 2 × iorder1ldsw<2iq+1+2×iorder-1,
orlwk < 2iq + 3lwk<2iq+3.

Accuracy

Not applicable.

Further Comments

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(ts))ds,  0tT,
0
y(t)=f(t)+1π0tK(t,s)y(s) t-s ds,  0tT,
(2)
using nag_inteq_abel_weak_weights (d05by). In (2), K(t,s)K(t,s) and f(t)f(t) are given and the solution y(t)y(t) is sought on a uniform mesh of size hh such that T = n × hT=n×h. Discretization of (2) yields
2p2 i
yi = f(i × h) + sqrt(h)Wi,jK(i × h,j × h)yj + sqrt(h)ωijK(i × h,j × h)yj,
j = 0 j = 2p1
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)yiy(i×h), for i = 1,2,,ni=1,2,,n. We propose the following algorithm for computing yiyi from (3) after a call to nag_inteq_abel_weak_weights (d05by):
(a) Set n = 2iq + 1 + 2 × iorder2n=2iq+1+2×iorder-2 and h = T / nh=T/n.
(b) Equation (3) requires 2 × iorder22×iorder-2 starting values, yjyj, for j = 1,2,,2 × iorder2j=1,2,,2×iorder-2, with y0 = f(0)y0=f(0). These starting values can be computed by solving the system
2 × iorder2
yi = f(i × h) + sqrt(h)sw(i + 1,j + 1)K(i × h,j × h)yj, i = 1,2,,2 × iorder2.
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 × iorder1,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,,ni=2×iorder-1,2×iorder,,n to compute yiyi from:
i1
(1sqrt( h )wt(1)K(i × h,i × h))yi = σi + sqrt(h)wt(ij + 1)K(i × h,j × h)yj.
j = 2 × iorder1
( 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



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