Symbol Error Probability of DF Relay Selection over Arbitrary Nakagami-𝑚 Fading Channels

Abstract

We present a new analytical expression for the moment generating function (MGF) of the end-to-end signal-to-noise ratio of dual-hop decode-and-forward (DF) relaying systems with relay selection when operating over Nakagami-𝑚 fading channels. The derived MGF expression, which is valid for arbitrary values of the fading parameters of both hops, is subsequently utilized to evaluate the average symbol error probability (ASEP) of 𝑀-ary phase shift keying modulation for the considered DF relaying scheme under various asymmetric fading conditions. It is shown that the MGF-based ASEP performance evaluation results are in excellent agreement with equivalent ones obtained by means of computer simulations, thus validating the correctness of the presented MGF expression.

1. Introduction

Cooperative communication through relay nodes has been shown to be capable of extending the radio coverage and improving the reliability of emerging wireless systems [1–3]. One of the bandwidth efficient dual-hop cooperative techniques combines the decode-and-forward (DF) relaying protocol with relay selection (RS) [4]. The performance of this technique has been studied in [4] over Rayleigh fading channels and in [5–7] for the more general Nakagami-𝑚 fading channel model. However, in the latter works analytical expressions for the moment generating function (MGF) of the end-to-end signal-to-noise ratio (SNR) have been presented, which are valid only for the special case where the Nakagami 𝑚-parameter of both hops takes integer values. In [5], the authors based their analysis on the tight approximation for the end-to-end SNR presented in [8], whereas [6, 7] utilized the end-to-end SNR characterization of [9]. Nevertheless, in realistic wireless communication scenarios, estimators for 𝑚 from field measurement data typically result in arbitrary 𝑚 values [10]. Moreover, restricting 𝑚 to only integer values severely limits the advantageous property of the Nakagami-𝑚 fading distribution to adequately approximate the Rice and Hoyt ones [11]. Very recently, based on [9], the authors in [12] investigated the error probability of opportunistic DF relaying over Nakagami-𝑚 fading channels with arbitrary 𝑚.

In this paper, capitalizing on the approach of [8], we present a new closed-form representation for the MGF of the end-to-end SNR of dual-hop DF RS-based systems which is valid for arbitrary-valued parameters for both Nakagami-𝑚 faded hops. In addition, the derived expression is utilized to evaluate the average symbol error probability (ASEP) of 𝑀-ary phase shift keying (PSK) modulation for the considered relaying scheme. Numerically evaluated ASEP results match perfectly with equivalent results obtained from computer simulations and clearly demonstrate that ASEP is rather sensitive to even slight variations of any of the hops’ fading parameter.

The remainder of this paper is organized as follows. Section 2 presents in brief the corresponding signal and system model. A new closed-form representation for the MGF of DF with RS over arbitrary Nakagami-𝑚 fading channels is derived in Section 3. An application of the derived expression in evaluating the ASEP performance of 𝑀-PSK is demonstrated in Section 4, while closing remarks are provided in Section 5.

2. Signal and System Model

Consider a wireless dual-hop DF cooperative system comprising a source node 𝑆, a destination node 𝐷, and 𝐿 half-duplex relays, each of which is denoted by 𝑅_𝑘, 𝑘=1,2,…,𝐿. During the first transmission time slot, 𝑆 broadcasts a symbol to the relays. In the second time slot, the relay with the most favorable channel to 𝐷 is selected to transmit, while 𝑆 remains silent. Node 𝐷 is assumed to possess perfect channel state information so that maximum-likelihood combining of the signals from 𝑆 and the selected relay can be employed.

The wireless channels between any pair of system nodes are assumed to be subject to the Nakagami-𝑚 fading distribution [13]. Furthermore, transmission is corrupted by additive white Gaussian noise (AWGN) with single-sided power spectral density 𝑁₀. The fading parameters for the 𝑆→𝐷, 𝑆→𝑅_𝑘, and 𝑅_𝑘→𝐷 channels are denoted by 𝑚₀, 𝑚_1𝑘, and 𝑚_2𝑘, respectively, whereas the corresponding instantaneous received SNRs as 𝛾₀, 𝛾_1𝑘, and 𝛾_2𝑘, respectively. It is noted that unlike many other studies [4–7], no restrictions are applied on the Nakagami-𝑚 fading parameters of both hops. In other words, 𝑚₀ and each 𝑚_ℓ𝑘, with ℓ=1,2 and 𝑘=1,2,…,𝐿, can take any arbitrary positive real value as long as it obeys the necessary condition for the Nakagami-𝑚 distribution, that is, 𝑚₀,𝑚_ℓ𝑘≥0.5. The average received SNRs at the first and second hop are defined as 𝛾₀=Ω_𝑆𝐷𝐸_𝑠/𝑁₀, 𝛾_1𝑘=Ω_𝑆𝑅𝑘𝐸_𝑠/𝑁₀, and 𝛾_2𝑘=Ω_𝑅𝑘𝐷𝐸_𝑠/𝑁₀, respectively, where Ω_𝐴𝐵 denotes the average fading power of the 𝐴→𝐵 channel and 𝐸_𝑠 is the energy of the symbol transmitted. Clearly, 𝛾₀ and 𝛾_ℓ𝑘’s are Gamma distributed random variables which are assumed to be statistically independent.

As shown in [8], a dual-hop transmission 𝑆→𝑅_𝑘→𝐷 employing the DF protocol can be modeled as a single-hop transmission 𝑆→𝐷 with the equivalent instantaneous received SNR at 𝐷 bounded as 𝛾_min−3.24𝑀<𝛾_𝑘≤𝛾_min. (1) In (1), 𝛾_min=min{𝛾_1𝑘,𝛾_2𝑘} and 𝑀 denotes the transmitted signal’s modulation order. For relatively large SNR and 𝑀 values, setting 𝛾_𝑘=𝛾_min yields a tight approximation for 𝛾_𝑘 [8]. To this effect, by using RS, the corresponding received SNR at 𝐷 can be accurately approximated as 𝛾_RS=max_{𝑘=1,2,…,𝐿}𝛾_𝑘 [5]. Utilizing next a time-diversity version of maximal-ratio combining [11], the instantaneous received SNR at 𝐷’s output can be expressed as 𝛾_𝐷=𝛾₀+𝛾_RS. (2)

3. MGF of DF RS in Nakagami-𝑚 Fading

By differentiating [5, equation (3)], a closed-form expression for the probability density function (PDF) of 𝛾_𝑘 for any arbitrary value of 𝑚_1𝑘,𝑚_2𝑘≥0.5 can be obtained as follows: 𝑓_𝛾𝑘(𝑥)=²_𝑝=1𝑎^𝑚𝑝𝑘_𝑝𝑘Γ𝑚_𝑝𝑘𝑥^{𝑚𝑝𝑘−1}exp−𝑎_𝑝𝑘𝑥×𝐺𝑚_{(3−𝑝)𝑘},𝑎_{(3−𝑝)𝑘}𝑥, (3) where 𝑎_1𝑘=𝑚_1𝑘/𝛾_1𝑘, 𝑎_2𝑘=𝑚_2𝑘/𝛾_2𝑘, and 𝐺(𝑐,𝑑𝑥)=Γ(𝑐,𝑑𝑥)/Γ(𝑐), with 𝑐,𝑑∈ℝ, while Γ(⋅) and Γ(⋅,⋅) denote the incomplete [14, equation (8.310)] and the upper incomplete Gamma function [14, equation (8.350.2)], respectively. Substituting [5, equation (3)] and (3) into [5, equation (6)] and after some algebraic manipulations, the following expression for the PDF of 𝛾_RS is deduced: 𝑓_𝛾RS(𝑥)=^2𝐿_𝑝=1^𝐿_𝑘=1𝑙=0∑^𝑙_𝑡=1(−1)^𝑙𝑎^𝑚𝑝𝑘_𝑝𝑘𝑙!Γ𝑚_𝑝𝑘𝑥^{𝑚𝑝𝑘−1}×exp−𝑎_𝑝𝑘𝑥𝐺𝑚_𝑝𝑛𝑡,𝑎_𝑝𝑛𝑡𝑥×𝐺𝑚_{(3−𝑝)𝑘},𝑎_{(3−𝑝)𝑘}𝑥𝐺𝑚_{(3−𝑝)𝑛𝑡},𝑎_{(3−𝑝)𝑛𝑡}𝑥, (4) where symbol ∑ is used for short-hand representation of multiple summations of the form ∑^𝐿_𝑛1=1∑^𝐿_𝑛2=1⋯∑^𝐿_𝑛𝑙=1 with 𝑛₁≠𝑛₂≠⋯≠𝑛_𝑙.

In order to derive an explicit expression for the MGF of 𝛾_R𝑆, which is defined as 𝑀_𝛾RS(𝑠)=∫^∞₀𝑓_𝛾RS(𝑥)exp(𝑠𝑥)𝑑𝑥 [11, equation (1.2)], for arbitrary valid values of both 𝑚_1𝑘 and 𝑚_2𝑘, one is required to analytically evaluate integrals that involve combinations of arbitrary powers, exponentials, and 𝐺(⋅,⋅) functions. To this end, by expressing all G(⋅,⋅)’s in (4) according to [14, equation (8.354.2)], utilizing [14, equation (8.310)], and performing some rather long but basic algebraic manipulations, one obtains the following explicit expression for 𝑀_𝛾RS(𝑠): 𝑀_𝛾RS(𝑠)=^𝐿𝐿_𝑘=1𝑙=0∑^𝑙_𝑡=1(−1)^𝑙𝑙!²_ℓ=11Γ𝑚_ℓ𝑘Γ𝑚_ℓ𝑛𝑡×𝑎^𝑚1𝑘8_1𝑘𝑝=1𝒜_𝑝+𝑎^𝑚2𝑘8_2𝑘𝑝=1ℬ_𝑝 (5) which is valid for 𝑎_1𝑘,𝑎_2𝑘≠𝑠 for all 𝑘=1,2,…,𝐿. In (5), parameters 𝒜_𝑝 and ℬ_𝑝, with 𝑝=1,2,…,8, are given by 𝒜₁ℬ₁=²_ℓ=1Γ𝑚_ℓ𝑛𝑡Γ(Ξ)𝐻₃(1,1,1),𝒜₂ℬ₂=−Γ𝑚_2𝑛𝑡Γ(Ξ)^∞_𝑖=0(−1)^𝑖𝑎𝑖!^{𝑚𝑡1𝑛+𝑖1𝑛𝑡}𝐻₂𝑚_1𝑛𝑡,𝒜+𝑖,1,1₃ℬ₃=−²_ℓ=1Γ𝑚_ℓ𝑛𝑡^∞_𝑖=0(−1)^𝑖𝑎𝑖!^2𝑘𝑎_1𝑘^Ξ+𝑖𝐻₂(),𝒜Ξ+𝑖,1,1₄ℬ₄=Γ𝑚_2𝑛𝑡^∞_𝑖,𝑗=0(−1)^𝑖+𝑗𝑎𝑖!𝑗!^{𝑚𝑡1𝑛+𝑗1𝑛𝑡}𝑎_2𝑘𝑎_1𝑘^Ξ+𝑖×𝐻₁Ξ+𝑖,𝑚_1𝑛𝑡,𝒜+𝑗,1₅ℬ₅=−Γ𝑚_1𝑛𝑡Γ(Ξ)^∞_𝑖=0(−1)^𝑖𝑎𝑖!^{𝑚𝑡2𝑛+𝑖2𝑛𝑡}𝐻₂𝑚_2𝑛𝑡,𝒜+𝑖,1,1₆ℬ₆=Γ(Ξ)^∞_𝑖,𝑗=0(−1)^𝑖+𝑗𝑎𝑖!𝑗!^{𝑚𝑡1𝑛+𝑖1𝑛𝑡}𝑎^{𝑚𝑡2𝑛+𝑗}_2𝑛𝑡×𝐻₁𝑚_1𝑛𝑡+𝑖,𝑚_2𝑛𝑡,𝒜+𝑗,1₇ℬ₇=Γ𝑚_1𝑛𝑡^∞_𝑖,𝑗=0(−1)^𝑖+𝑗𝑎𝑖!𝑗!^2𝑘𝑎_1𝑘^Ξ+𝑖𝑎^{𝑚𝑡2𝑛+𝑗}_2𝑛𝑡×𝐻₁Ξ+𝑖,𝑚_2𝑛𝑡,𝒜+𝑗,1₈ℬ₈=−^∞_{𝑖,𝑗,𝑧=0}(−1)^{𝑖+𝑗+𝑧}𝑎𝑖!𝑗!𝑧!^2𝑘𝑎_1𝑘^Ξ+𝑖𝑎^{𝑚𝑡1𝑛+𝑗}_1𝑛𝑡𝑎^{𝑚𝑡2𝑛+𝑧}_2𝑛𝑡×𝐻₀Ξ+𝑖,𝑚_1𝑛𝑡+𝑗,𝑚_2𝑛𝑡+𝑧. (6)(7) In (6)-(7) the notation ^𝑥_𝑦={^𝑐_𝑑} is used for short-hand representation of 𝑥=𝑐 and 𝑦=𝑑, and symbols Ξ and Υ represent 𝑚_2𝑘𝑚_1𝑘 and 𝑎_1𝑘𝑎_2𝑘, respectively, whereas the function of 𝑠 denoted by 𝐻_𝑛(𝜅,𝜆,𝜇), with 𝜅,𝜆,𝜇∈ℝ and 𝑛=0, 1, 2 and 3, is defined as 𝐻_𝑛=Γ𝜅,𝜆,𝜇^{𝑚1𝑘𝑚2𝑘}+𝜅+𝜆+𝜇−𝑛𝜅𝜆𝜇(Υ−𝑠)^{𝑚1𝑘𝑚2𝑘+𝜅+𝜆+𝜇−𝑛}. (8)

Importantly, after careful inspection of 𝒜_𝑝’s and ℬ_𝑝’s for 𝑝=2,3,…,8 in (5), it follows that they can all be expressed in terms of well-known generalized hypergeometric functions [14, equation (9.14)]. To this end, for the 𝒜_𝑝’s and ℬ_𝑝’s with 𝑝=2, 3, and 5, by expressing every term of the form (𝑥+𝑦) as (𝑥+𝑦)!/(𝑥+𝑦−1)! and using the identity Γ(𝑥+𝑦)=(𝑥)_𝑦Γ(𝑦), where (𝑥)_𝑦 is the Pochhammer symbol [14, page xliii], yields 𝒜₂ℬ₂=Γ𝑚_2𝑛𝑡Γ(Ξ)𝐼𝑚_1𝑛𝑡,𝑎_1𝑛𝑡,𝒜₃ℬ₃=²_ℓ=1Γ𝑚_ℓ𝑛𝑡𝐼𝑎Ξ,^2𝑘𝑎_1𝑘,𝒜₅ℬ₅=Γ𝑚_1𝑛𝑡Γ(Ξ)𝐼𝑚_2𝑛𝑡,𝑎_2𝑛𝑡, (9) where 𝐼(𝜅,𝜆) is given by 𝐼(𝜅,𝜆)=−𝜆^𝜅𝐻₂(𝜅,1,1)₂𝐹₁𝑚_1𝑘𝑚_2𝑘𝜆+𝜅,𝜅;𝜅+1;.𝑠−Υ (10) In (10), ₂𝐹₁(⋅) represents the Gauss hypergeometric function [14, equation (9.14.2)]. Likewise, for the 𝒜_𝑝’s and ℬ_𝑝’s with 𝑝=4, 6, and 7, each infinite series can be ultimately expressed in terms of the Appell hypergeometric function 𝐹₂(⋅) [14, equation (9.180.2)] as 𝒜₄ℬ₄=Γ𝑚_2𝑛𝑡𝐽Ξ,𝑚_1𝑛𝑡,𝑎_2𝑘𝑎_1𝑘,𝑎_1𝑛𝑡,𝒜₆ℬ₆=Γ(Ξ)𝐽𝑚_1𝑛𝑡,𝑚_2𝑛𝑡,𝑎_1𝑛𝑡,𝑎_2𝑛𝑡,𝒜₇ℬ₇=Γ𝑚_1𝑛𝑡𝐽Ξ,𝑚_2𝑛𝑡,𝑎_2𝑘𝑎_1𝑘,𝑎_2𝑛𝑡, (11) where 𝐽(𝜅,𝜆,𝜇,𝜈), with 𝜈∈ℝ, is given by 𝐽𝜅,𝜆,𝜇,𝜈=𝜇^𝜅𝜈^𝜆𝐻₁(𝜅,𝜆,1)×𝐹₂𝑚_1𝑘𝑚_2𝑘𝜇+𝜅+𝜆;𝜅,𝜆;𝜅+1,𝜆+1;,𝜈𝑠−Υ.𝑠−Υ (12) Finally, by using once more all aforementioned identities, it follows that the 𝒜₈ and ℬ₈ coefficients can be expressed in terms of the generalized Lauricella function 𝐹_𝐴⁽³⁾(⋅) [15, equation (1.1)], yielding 𝒜₈ℬ₈=−𝑎^𝜆_1𝑛𝑡𝑎^𝜇_2𝑛𝑡𝑎_2𝑘𝑎_1𝑘^Ξ𝐻₀Ξ,𝜆,𝜇×𝐹_𝐴^{(3)𝑚1𝑘𝑚2𝑘+Ξ+𝜆+𝜇;Ξ,𝜆,𝜇;}_{Ξ+1,𝜆+1,𝜇+1;𝑎2𝑘𝑎1𝑘}𝑠−Υ^,𝑎_1𝑛𝑡𝑠−Υ^,𝑎_2𝑛𝑡.𝑠−Υ (13) To this effect, substituting (6), (9), (11), and (13) to (5), a novel expression for 𝑀_𝛾RS(𝑠) is deduced which is valid for arbitrary 𝑚_1𝑘,𝑚_2𝑘≥0.5. It is noted that in Appendices A and B we present MATLAB routines for computational efficient implementations of functions 𝐹₂(⋅) and 𝐹³_𝐴(⋅), respectively.

Using the previously derived 𝑀_𝛾RS(𝑠) expression, the MGF of 𝛾₀ and by recalling that 𝛾_RS and 𝛾₀ are statistically independent, the MGF of 𝛾_𝐷 is straightforwardly deduced, namely, 𝑀_𝛾𝐷(𝑠)=𝑚₀𝛾₀^𝑚0𝑚𝑠+⁰𝛾₀^−𝑚0𝑀_𝛾RS(𝑠). (14)

4. ASEP Performance Evaluation

The ASEP performance of various modulation schemes for the considered dual-hop DF relaying system with RS over Nakagami-𝑚 fading channels can be directly evaluated using the 𝑀_𝛾𝐷(𝑠) expression and the MGF-based approach presented in [11, Chapter 1]. For example, the ASEP of 𝑀-PSK modulation is easily obtained as 𝑃_se=1𝜋₀^{𝜋−𝜋/𝑀}𝑀_𝛾𝐷sin²(𝜋/𝑀)sin²(𝜃)𝑑𝜃. (15)

Let us assume a relaying system with 𝐿=3 relays, no direct link 𝑆→𝐷, and the following four asymmetric fading scenarios: (i) Scenario A: 𝑚₁₁=Ω₁₁=0.6, 𝑚₁₂=Ω₁₂=1.1, 𝑚₁₃=Ω₁₃=1.7, and 𝑚_2𝑘=Ω_2𝑘=0.7 for 𝑘=1, 2, and 3; (ii) Scenario B: 𝑚₁₁=Ω₁₁=1.1, 𝑚₁₂=Ω₁₂=2.1, 𝑚₁₃=Ω₁₃=3.1, and 𝑚_2𝑘=Ω_2𝑘=1.1 for 𝑘=1, 2, and 3; (iii) Scenario C: 𝑚₁₁=Ω₁₁=1.1, 𝑚₁₂=Ω₁₂=2.3, 𝑚₁₃=Ω₁₃=3.8, and 𝑚_2𝑘=Ω_2𝑘=1.8 for 𝑘=1, 2, and 3; and (iv) Scenario D: 𝑚₁₁=Ω₁₁=1.6, 𝑚₁₂=Ω₁₂=3.3, 𝑚₁₃=Ω₁₃=4.2, and 𝑚_2𝑘=Ω_2𝑘=2.6 for 𝑘=1, 2 and 3. The ASEP performance, 𝑃_se, of 8-, 16- and 32-PSK modulation is depicted in Figures 1, 2, and 3, respectively, as a function of the average symbol to noise power, 𝐸_𝑠/𝑁₀, over various Nakagami-𝑚 fading conditions. As expected, 𝑃_se improves with increasing 𝐸_𝑠/𝑁₀ and/or decreasing 𝑀 and/or increasing any of the fading parameters. Furthermore, all figures clearly depict the excellent agreement between the numerically evaluated ASEP results and the equivalent ones obtained from Monte Carlo simulations. In addition, it is observed that the higher is the 𝐸_𝑠/𝑁₀, the more sensitive is 𝑃_se to slight variations of the fading conditions. More specifically, for the case of 8-PSK modulation, it can be observed that the differences in the 𝑃_se curves among the four considered fading scenarios for the low 𝐸_𝑠/𝑁₀ regime, for example, at 5dB, are 51% between Scenarios A and B; 52% between Scenarios B and C; and 54% between Scenarios C and D. For the same modulation order, the differences in the high 𝐸_s/𝑁₀ regime, for example, at 20dB, are: 98% between Scenario A and B; 100% between Scenario B and C; and 103% between Scenario C and D. Likewise, for the case of 16-PSK and for 𝐸_𝑠/𝑁₀=5dB, the differences in the 𝑃_se curves for the different scenarios are 23% between Scenarios A and B; 29% between Scenarios B and C; and 36% between Scenarios C and D. For 16-PSK modulation and for 𝐸_𝑠/𝑁₀=20dB, the differences among the various 𝑃_se curves are 84% between Scenarios A and B; 90% between Scenarios B and C; and 98% between Scenarios C and D. Finally, for 32-PSK modulation, it is evident that the 𝑃_se curves differ at 𝐸_𝑠/𝑁₀=5dB: 11% between Scenarios A and B; 14% between Scenarios B and C; 19% between Scenarios C and D whereas, for 𝐸_𝑠/𝑁₀=20dB, the 𝑃_se differences among scenarios are: 46% between Scenarios A and B; 62% between Scenarios B and C; and 79% between Scenario C and D. It is evident from the above quantitative results that the sensitivity of the ASEP on the fading severity parameter 𝑚 is high in all considered communication scenarios and particularly in the case of modulation schemes with low 𝑀 and/or high 𝐸_𝑠/𝑁₀ values.

Figure 1 

Average symbol error probability, 𝑃se, of 8-PSK of dual-hop DF relaying with RS as a function of the average symbol to noise power, 𝐸𝑠/𝑁0, for 𝐿=3 relays and various Nakagami-𝑚 fading conditions.

Figure 2 

Average symbol error probability, 𝑃se, of 16-PSK of dual-hop DF relaying with RS as a function of the average symbol to noise power, 𝐸𝑠/𝑁0, for 𝐿=3 relays and various Nakagami-𝑚 fading conditions.

Figure 3 

Average symbol error probability, 𝑃se, of 32-PSK of dual-hop DF relaying with RS as a function of the average symbol to noise power, 𝐸𝑠/𝑁0, for 𝐿=3 relays and various Nakagami-𝑚 fading conditions.

5. Conclusion

This work was devoted to the derivation of a new analytical expression for the MGF of the end-to-end SNR of dual-hop DF relaying communication systems with RS over Nakagami-𝑚 fading conditions. The presented expression involves well-known generalized hypergeometric functions and is valid for any arbitrary value of the fading parameters of both hops. Using the MGF-based approach, the ASEP of 𝑀-PSK modulation for the considered system was evaluated and a perfect match with equivalent computer simulated performance results was shown. More importantly, it was evident that the ASEP is sensitive to even slight variations of any of the hops’ fading parameters and particularly for low order modulation schemes and high SNR values.

Appendices

A. MATLAB Program for the Evaluation of 𝐹₂(⋅)

The Appell hypergeometric function 𝐹₂(⋅) is defined by the double infinite series given by [14, equation (9.180.2)]. An equivalent expression for 𝐹₂(⋅) in a double integral form with finite limits is given by [14, equation (9.184.2)]. The latter expression can be easily evaluated using the MATLAB function shown in Algorithm 1.

Function  F2 = Appell(a,b1,b2,c1,c2,x,y); f1 = gamma(c1).*gamma(c2); f2 = gamma(b1).*gamma(b2); f3 = gamma(c1-b1).*gamma(c2-b2); f = f1./(f2.*f3); Q = @(u,v)f.*u.̂‌(b1-1).*(v.̂‌(b2-1)).*...    ((1-u).̂(c1-b1-1)).*...    ((1-v).̂(c2-b2-1))./...    ((1-u.*x-v.*y).̂a); F2 = dblquad(Q,0,1,0,1);

Algorithm 1 

MATLAB Program for evaluating 𝐹₂(⋅).

B. MATLAB Program for the Evaluation of 𝐹_𝐴⁽³⁾(⋅)

The generalized Lauricella function 𝐹_𝐴⁽³⁾(⋅) is expressed according to the triple infinite series given by [15, equation (1.1)]. With the aid of the integral representation for 𝐹₂(⋅) in [14, equation (9.184.2)], the following triple integral form for 𝐹_𝐴⁽³⁾(⋅) with finite limits is deduced: 𝐹_𝐴⁽³⁾𝑎;𝑏₁,𝑏₂,𝑏₃;𝑐₁,𝑐₂,𝑐₃=;𝑥,𝑦,𝑧₁₀₁₀₁₀𝑓(𝑢,𝑣,𝑤)𝑑𝑢𝑑𝑣𝑑𝑤, (B.1) where 𝑎,𝑏₁,𝑏₂,𝑏₃,𝑐₁,𝑐₂,𝑐₃,𝑥,𝑦,𝑧∈ℝ and 𝑓()=𝑢𝑢,𝑣,𝑤^𝑏1−1𝑣^𝑏2−1𝑤^𝑏3−1(1−𝑢)^{𝑐1𝑏1−1}Γ𝑏₁Γ𝑏₂Γ𝑏₃𝑐_1−𝑏1𝑐_2−𝑏2×(1−𝑣)^{𝑐2−𝑏2−1}(1−𝑤)^{𝑐3−𝑏3−1}𝑐_3−𝑏31−𝑢𝑥−𝑣𝑦−𝑤𝑧^𝑎. (B.2) To this effect, the MATLAB function shown in Algorithm 2 evaluates 𝐹_𝐴⁽³⁾(⋅) using (B.1).

Function  FA = Lauricella  (a,b1,b2,b3,c1,c2,c3,x,y,z); f1 = gamma(c1).*gamma(c2).*gamma(c3); f2 = gamma(b1).*gamma(b2).*gamma(b3); f3 = gamma(c1-b1).*gamma(c2-b2).*gamma(c3-b3); f = f1./(f2.*f3); Q = @(u,v,w)f.*(u.̂‌(b1-1)).*...    (v.̂‌(b2-1)).*(w.̂‌(b3-1)).*...    ((1-u).̂‌(c1-b1-1)).*...    ((1-v).̂‌(c2-b2-1)).*...    ((1-w).̂‌(c3-b3-1))./...    ((1-u.*x-v.*y-w.*z).̂‌a); FA = triplequad(Q,0,1,0,1,0,1);

Algorithm 2 

MATLAB Program for evaluating 𝐹_𝐴⁽³⁾(⋅).

Acknowledgments

The work of Dr. George C. Alexandropoulos has been partly supported by the General Secretarial Research and Technology (G.S.R.T.), Hellas, under the ARTEMIS research program R3-COP with the Project number 100233. The work of Dr. Khuong Ho-Van is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant number 102.04-2012.39.