11458 FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 36, No 3, September 2023, pp. 395-410 https://doi.org/10.2298/FUEE2303395S © 2023 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper HAAR WAVELET OPERATIONAL MATRIX BASED NUMERICAL INVERSION OF LAPLACE TRANSFORM FOR IRRATIONAL AND TRANSCENDENTAL TRANSFER FUNCTIONS Zdravko Stanimirović, Ivanka Stanimirović, Slobodanka Galović, Katarina Đorđević, Edin Suljovrujić Vinča Institute of Nuclear Sciences, National Institute of the Republic of Serbia, University of Belgrade Abstract. Irrational and transcendental functions can often be seen in signal processing or physical phenomena analysis as consequences of fractional-order and distributed- order models that result in fractional or partial differential equations. In cases when finding solution in analytical form tends to be difficult or impossible, numerical calculations such as Haar wavelet operational matrix method can be used. Haar wavelet establishes a direct procedure for transfer function inversion using the wavelet operational matrix for orthogonal function set integration. In this paper an inverse Laplace transform of irrational and transcendental transfer functions using Haar wavelet operational matrix is proposed. Results for a number inverse Laplace transforms are numerically solved and compared with the analytical solutions and solutions provided by commonly used Invlap and NILT algorithms. This approach is useful when the original cannot be represented by an analytical formula and validity of the obtained result needs to be crosschecked and error estimated. Key words: Haar wavelet, Laplace transform, maximum resolution level, numerical inversion, operational matrix 1. INTRODUCTION Laplace transform is a powerful tool in solving various problems in engineering and science. It usually simplifies differential equations that come from these areas. The inverse Laplace transformation is usually done with the use of Laplace transform tables in combination with simple algebraic manipulations. However, there are numerous cases when finding solution in analytical form tends to be difficult or impossible [1-6]. Irrational and transcendental functions can often be seen in signal processing or physical phenomena Received January 11, 2023; revised Fenruary 28, 2023; accepted March 13, 2023 Corresponding author: Zdravko Stanimirović Vinča Institute of Nuclear Sciences, National Institute of the Republic of Serbia, University of Belgrade, Serbia E-mail: zdravko.stanimirovic@vinca.rs 396 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ analysis as consequences of fractional-order and distributed-order models that result in fractional or partial differential equations [7-10]. Over the years, many numerical inverse Laplace transform algorithms are being used depending on the field of interest. Algorithms like Direct, Invlap, Weeks, Gavsteh, NILT, etc. differ in parameter choice, effectiveness, computational speed and reliability [11]. Lately, researchers are often using Haar wavelet operational matrix for finding inverse Laplace transforms [12-14]. Haar wavelet establishes a direct procedure for transfer function inversion using the wavelet operational matrix for orthogonal function set integration. This approach is straightforward and suitable for computer programming using special purpose programming environments that allow matrix manipulations. However, when irrational and transcendental functions are in question problems may occur in case of fractional power of the matrix, matrices that appear as arguments of transcendental functions, etc. The objective of this paper is to evaluate the application potential of Haar wavelet operational matrix for finding irrational and transcendental inverse Laplace transforms. The effectiveness of the method is estimated by comparing obtained results with analytical ones and results obtained by NILT and Invlap algorithms. Influence of the maximum resolution level of the operational matrix on the agreement with the analytical solution is examined. Also, standard and absolute errors are calculated and analyzed as well as the orders of convergence at points of sharp turns. 2. HAAR WAVELET METHOD Haar functions are defined in the interval of [0,τ) [12-13] by: ℎ0(𝑡) = 𝑚− 1 2 and ℎ1(𝑡) = 𝑚− 1 2 { 1, 2−𝑗𝜏(𝑘 − 1) ≤ 𝑡 < 2−𝑗𝜏(𝑘 − 1 2⁄ ) −1, 2−𝑗𝜏(𝑘 − 1 2⁄ ) ≤ 𝑡 < 2−𝑗𝜏𝑘 0, elsewere in [0, τ) (1) where i= 0,1,2…, (m-1), m=2α, αϵZ+, m being denoted as the maximum level of resolution. Integer decomposition of the index i is designated by j and k (e.g., i=2j+k-1 where k=1,2,…, 2j). Scaling function h0(t) is a constant, while the fundamental square wave h1(t) is known as the mother wavelet. All other wavelets are generated from the fundamental one: ℎ𝑖(𝑡) = 𝑚− 1 2 { 2𝑗/2, 2−𝑗𝜏(𝑘 − 1) ≤ 𝑡 < 2−𝑗𝜏(𝑘 − 1 2⁄ ) −2𝑗/2, 2−𝑗𝜏(𝑘 − 1 2⁄ ) ≤ 𝑡 < 2−𝑗𝜏𝑘 0, otherwise in [0, τ) (2) Any function x(t) can be expanded into a Haar series: 𝑥(𝑡) = ∑ 𝑐𝑖 𝑚−1 0 ℎ𝑖(𝑡), 0 ≤ 𝑡 < τ (3) with Haar coefficient ci: 𝑐𝑖 = 𝑚 ∫ 𝑥(𝑡)ℎ𝑖(𝑡)𝑑𝑡. 𝑡 0 (4) In the matrix form equation (3) becomes: 𝑥𝑇 = 𝑐𝑇 ∙ 𝐻 (5) where x designates the discrete form of the function x(t). cT is the Haar coefficient vector: 𝑐𝑇 = [𝑐0 𝑐1 … 𝑐𝑚−1] (6) Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 397 Haar function vector is defined as: 𝐻 = [ℎ0 ℎ1 ⋯ ℎ𝑚−1]𝑇 (7) Integration of the Haar wavelet function H can be written as: ∫ 𝐻(𝑡)𝑑𝑡 = 𝑄𝐻 ∙ 𝐻 𝑖 0 , 0 ≤ 𝑡 < 𝜏 (8) where QH is the Haar operational matrix for integration. If we take the H(t)=H·Bm(t) form of Haar wavelet function where Bm(t) is the block pulse function, then: ∫ 𝐻(𝑡)𝑑𝑡 = ∫ 𝐻 ∙ 𝐵𝑚(𝑡)𝑑𝑡 = 𝑖 0 𝐻 ∙ 𝑄𝑚 ∙ 𝐵𝑚(𝑡) 𝑖 0 (9) The operational matrix for integration of the block pulse function Qm in the interval of [0,1) is given by: 𝑄𝑚 = (2 ∙ 𝑚)−1 [ 1 2 ⋯ 2 0 1 ⋯ ⋮ ⋮ 0 0 ⋯ ⋱ 2 0 1 ] 𝑚×𝑚 , 0 ≤ 𝑡 < 1 (10) and in the interval of [0,τ): 𝑄𝑚 = (2 ∙ 𝑚)−1 [ 2 ∙ 𝑚 ∙ 𝑄𝑚/2 −𝜏 ∙ 𝐻𝑚/2 𝑇 𝜏 ∙ 𝐻𝑚/2 𝑇 0𝑚/2 ] , 0 ≤ 𝑡 < 𝜏 (11) Since Bm(t) can be neglected as an identity matrix, the generalized Haar operational matrix can be obtained from (8) and (9): 𝑄𝐻 = 𝐻 ∙ 𝑄𝑚 ∙ 𝐻𝑇 (12) 3. NUMERICAL INVERSION OF LAPLACE TRANSFORM USING THE OPERATIONAL MATRIX In order to determine inverse Laplace transform of the function x(t) using the operational matrix of integration we should start from the fact that integration in time domain (∫ 𝑥(𝜏)𝑑𝜏 𝑖 0 ) corresponds to the multiplication of 1/s in the s domain ( 1 𝑠 ∙ 𝑋(𝑠)). The following transfer function: 𝑋(𝑠) = 𝑏 𝑠+𝑎 = 𝑏 𝑠 1+ 𝑎 𝑠 = �̂� ( 1 𝑠 ) (13) is a solution to the following equation given in a discrete form: 𝑏 ∙ [1 1 … 1] = 𝒙 + 𝑎 ∙ ∫ 𝒙𝑑𝑡 𝑖 0 (14) where the inversion of the Laplace transform column vector xT can be replaced by a product of coefficient vector cT and Haar wavelet matrix H: 𝑏 ∙ [1 1 … 1] ∙ 𝐻−1 =𝑐𝑇 ∙ [𝐼 + 𝑎 ∙ 𝑄𝐻] (15) In eq. (14) I is the m×m identity matrix and the matrix of integration QH replaced the integral sign. Based on the eq. (12) the discrete form of the transfer function X(s) can be expressed as: 398 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ �̂�(𝑄𝐻) = [𝑏 ∙ 𝑄𝐻] ∙ [𝐼 + 𝑎 ∙ 𝑄𝐻]−1 (16) From (9) and (14) the coefficient vector becomes: 𝑐𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ �̂�(𝑄𝐻) (17) Accordingly, the inversion of the Laplace transform X(s) can be given by: 𝑥𝑇 = 𝑐𝑇 ∙ 𝐻 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ �̂�(𝑄𝐻) ∙ 𝐻 (18) 4. NUMERICAL INVERSE LAPLACE TRANSFORM ALGORITHMS Over the years several numerical inverse Laplace transform algorithms based on different numerical methods have been developed. Weeks algorithm is based on bilinear transformations and the Laguerre polynomial series expansion [15]. Another algorithm for inverting the Laplace transform is Direct. Direct is based on the application of the trapezoidal rule to the Bromwich integral [16]. Invlap method is based on De Hoog’s algorithm [17]. Gavsteh is a numerical inverse Laplace transform algorithm based on Gaver-Stehfest method [18]. Fast Fourier transformation is the base of the NILT algorithm. In combination with the ε-algorithm it provides high speed computation and accuracy [19]. There is also NILT based on the Zakian method [20] - method that is based on the Fourier series method with Padé approximation. Application of numerical inverse Laplace transform algorithms was previously evaluated [11, 16, 19] and obtained results showed that a single method cannot provide optimum results for all purposes. For that reason, effectiveness of the Haar wavelet operational matrix for finding inverse Laplace transforms will be evaluated by comparison of results obtained by Haar wavelet method, analytical results and results obtained by NILT and Invlap algorithms. Invlap and Zakian method based NILT numerical inverse Laplace transform algorithms were chosen because they are effective and can deal with irrational and transcendental functions. This approach is especially useful when the analytical formula is not available and inverse Laplace transform can be found only numerically. In that case validity of the obtained result can be crosschecked and error can be estimated. 5. NUMERICAL INVERSION OF IRRATIONAL AND TRANSCENDENTAL LAPLACE TRANSFORM 5.1. Function 𝑿(𝒔) = 𝒃𝟐−𝒂𝟐 (𝒔−𝒂𝟐)(√𝒔+𝒃) In terms of 1/s transfer function becomes: �̂� ( 1 𝑠 ) = 1 𝑠√𝑠 𝑏2−𝑎2 (1− 𝑎2 𝑠 )(1+ 𝑏 √𝑠 ) (19) Each 1/s is then replaced by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = (𝑏2 − 𝑎2)𝑄𝐻 1.5[(𝐼 − 𝑎2𝑄𝐻) ∙ (𝐼 + 𝑏𝑄𝐻 0.5)]−1 (20) Then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ (𝑏2 − 𝑎2) ∙ 𝑄𝐻 1.5 ∙ [(𝐼 − 𝑎2𝑄𝐻) ∙ (𝐼 + 𝑏𝑄𝐻 0.5)]−1 ∙ 𝐻 (21) Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 399 In the code, fractional power of matrix is calculated indirectly using principal matrix logarithm where the matrix function is built on the principal scalar logarithm. The analytical inverse Laplace transform of the transfer function given in the heading 5.1 is: 𝑥(𝑡) = 𝑒𝑎2𝑡[𝑏 − 𝑎 𝑒𝑟𝑓(𝑎√𝑡)] − 𝑏𝑒𝑏2𝑡 𝑒𝑟𝑓𝑐(𝑏√𝑡) (22) Error function erf(x) and complementary error function erfc(x) are two the most widely used functions in science. These functions occur extensively in problems relating to heat conduction, diffusion and probability. In case of a=0.5 and b=1, for transfer function X(s), analytical result and numerical results obtained by Haar wavelet method (m=1024), Invlap and NILT algorithms are shown in Fig. 1(a) for interval [0,1) and in Fig. 1(c) for expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels (m=64, 256 and 1024) as well as for Invlap and NILT algorithms are presented in Fig. 1(b) for interval [0,1) and in Fig. 1(d) for expanded interval [0,τ). (a) (b) (c) (d) Fig. 1 Transfer function 𝑋(𝑠) = 𝑏2−𝑎2 (𝑠−𝑎2)(√𝑠+𝑏) inverse Laplace transform obtained analytically and numerically by Haar wavelet method (m=1024), Invlap and NILT algorithms for (a) interval [0,1) and (c) expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels and Invlap and NILT algorithms for intervals (b) [0,1) and (d) [0,τ) 400 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ For the transfer function given in the heading 5.1, in both [0,1) and [0,τ) interval, Haar wavelet method shows a good agreement with analytical solution as well as Invlap and NILT. When Haar wavelet method is in question, absolute error decreases with time during the [0,1) interval with minimum of the order of 10-7 for maximum m value. Standard error is of the order 10-1 for all three maximum resolution levels. When expanded interval is in question, after the [0,1) interval, Haar absolute error values fluctuate and have minimum values around t=3 when the increase in error with time starts. Standard error reaches the order of 101 while Invlap and NILT perform slightly better with standard errors of order of 100. Absolute error of NILT algorithm is almost constant during the whole time and is of order of 10-2. Absolute error of Invlap increases with time. It is in 10-10 – 10-1 range for [0,25) interval. 5.2. Function 𝑿(𝒔) = 𝟏 √𝒔∙(√𝒔+𝒂) In terms of 1/s transfer function becomes: �̂� ( 1 𝑠 ) = 1 𝑠 ∙ 1 (1+ 𝑎 √𝑠 ) (23) Each 1/s is then replaced by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 𝑄𝐻 ∙ (𝐼 + 𝑄𝐻 0.5 ∙ 𝑎)−1 (24) Then, the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 𝑄𝐻 ∙ (𝐼 + 𝑄𝐻 0.5 ∙ 𝑎)−1 ∙ 𝐻 (25) The analytical inverse Laplace transform of the transfer function given in the heading 5.2 is: 𝑥(𝑡) = 𝑒𝑎2𝑡 ∙ 𝑒𝑟𝑓𝑐(𝑎√𝑡) (26) In case of a=1, for transfer function X(s), analytical result and numerical results obtained by Haar wavelet method (m=1024), Invlap and NILT algorithms are shown in Fig. 2(a) for interval [0,1) and in Fig. 2(c) for expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels (m=64, 256 and 1024) as well as for Invlap and NILT algorithms are presented in Fig. 2(b) for interval [0,1) and in Fig. 2(d) for expanded interval [0,τ). For the transfer function given in the heading 5.2, for both intervals, Haar wavelet method shows a good agreement with analytical solution as well as Invlap and NILT. When Haar wavelet method is in question, absolute error decreases with time over the entire time span. It is in 10-2 – 10-6 range for maximum m. Standard error is of the order 10-1 for all three maximum resolution levels during the whole interval while Invlap and NILT perform better with standard errors of order of 10-3. Absolute error of NILT and Invlap algorithms are almost constant during the whole time. Haar wavelet method performs better than NILT whose absolute error is of the order of 10-2. Absolute error of Invlap algorithm is of the order of 10-9. Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 401 5.3. Function 𝑿(𝒔) = 𝒂𝒓𝒄𝒕𝒂𝒏 ( 𝒌 𝒔 ) When 1/s is replaced by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 𝑎𝑡𝑎𝑛(𝑘 ∙ 𝑄𝐻) (27) Then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 𝑎𝑡𝑎𝑛(𝑘 ∙ 𝑄𝐻) ∙ 𝐻 (28) Because of the matrix calculus, tangent function in the code was expressed using complex logarithms. The analytical inverse Laplace transform of the transfer function given in the heading 5.3 is: 𝑥(𝑡) = 1 𝑡 𝑠𝑖𝑛 (𝑘 ∙ 𝑡) (29) (a) (b) (c) (d) Fig. 2 Тransfer function 𝑋(𝑠) = 1 √𝑠∙(√𝑠+𝑎) inverse Laplace transform obtained analytically and numerically by Haar wavelet method (m=1024), Invlap and NILT algorithms for (a) interval [0,1) and (c) expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels and Invlap and NILT algorithms for intervals (b) [0,1) and (d) [0,τ). 402 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ This example is chosen because sinusoidal functions appear everywhere, and they play an important role in circuit analysis. Apart from electrical engineering they are seen in various branches of science and engineering. In case of k=1, for transfer function X(s), analytical result and numerical results obtained by Haar wavelet method (m=1024), Invlap and NILT algorithms are shown in Fig. 3(a) for interval [0,1) and in Fig. 3(c) for expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels (m=64, 256 and 1024) as well as for Invlap and NILT algorithms are presented in Fig. 3(b) for interval [0,1) and in Fig. 3(d) for expanded interval [0,τ). For the transfer function given in the heading 5.3, for both intervals, Haar wavelet method shows a good agreement with analytical solution as well as Invlap algorithm. When Haar wavelet method is in question, absolute error is approximately of the order of 10-7 during the [0,1) interval for maximum m. Haar wavelet method performs better than NILT whose absolute error is of the order of 10-5. Absolute error of Invlap algorithm is of the order of 10-9 for the same interval. During the expanded interval Haar wavelet absolute error is of the order of 10-4. Absolute error of Invlap algorithm is of the order of 10-10, while absolute error of NILT algorithm is of the order of 10-5 and increases to 10-1 at (a) (b) (c) (d) Fig. 3 Тransfer function 𝑋(𝑠) = 𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝑘 𝑠 ) inverse Laplace transform obtained analytically and numerically by Haar wavelet method (m=1024), Invlap and NILT algorithms for (a) interval [0,1) and (c) expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels and Invlap and NILT algorithms for intervals (b) [0,1) and (d) [0,τ). Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 403 t=11. Haar wavelet method standard error is of the order 10-2 for interval [0,1) and 10-1 for expanded interval [0,τ) for all three maximum resolution levels. Standard errors of Invlap and NILT are of order of 10-3 over the entire time span. 5.4. Function 𝑿(𝒔) = 𝟏 𝒔√𝒔 (𝟏 − 𝒆−𝑻𝒔) When we replace 1/s by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 𝑄𝐻 1.5 ∙ (𝐼 − 𝑒−𝑇𝑄𝐻 −1 ) (30) then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 𝑄𝐻 1.5 ∙ (𝐼 − 𝑒−𝑇𝑄𝐻 −1 ) ∙ 𝐻 (31) The analytical inverse Laplace transform of the transfer function given in the heading 5.4 is: 𝑥(𝑡) = 2 √𝜋 ∙ √𝑥 ∙ 𝐻(𝑥) − √𝑥 − 𝑇 ∙ 𝐻(𝑥 − 𝑇) (32) The Heaviside unit step function is used in the signal processing. It represents signals that switch on at specified times and stay switched on indefinitely. It is also used in structural mechanics to describe different structural loads, in engineering where periodic functions are represented, in physics for sudden changes (when breaks are being applied or during collisions), etc. In case of T=0.5 for transfer function X(s), analytical result and numerical results obtained by Haar wavelet method (m=1024), Invlap and NILT algorithms are shown in Fig. 4(a) for interval [0,1) and in Fig. 4(c) for expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels (m=64, 256 and 1024) as well as for Invlap and NILT algorithms are presented in Fig. 4(b) for interval [0,1) and in Fig. 4(d) for expanded interval [0,τ). When Haar wavelet method is in question, functions with sharp turns are challenging. For the transfer function given in the heading 5.4, over the [0,1) interval Haar method performs well at the sharp turn, almost as well as Invlap algorithm. Absolute errors at the peak for all algorithms are of the order of 10-2. However, because of the same maximum resolution level value (1024) and a longer period of time, Haar wavelet has poorer performances at the sharp turn than both Invlap and NILT in the [0,τ) interval. Over the [0,τ) interval absolute error of Harr wavelet method varied in the 10-2-10-4 range. For NILT algorithm, the absolute error has the highest value. It is of the order of 10-2. Apart from the sharp turn, absolute error of Invlap algorithm is of the order of 10-10. Haar wavelet method standard error is of the order 10-1 during the [0,1) interval and 10-2 during the [0,τ) interval for all three maximum resolution levels. Standard errors of Invlap and NILT are of order of 10-3 over the whole time span. From Figs. 1-4 can be seen that numerical inversion Laplace transform using Haar wavelet operational matrix performs very well in case of irrational and transcendental functions. Results obtained by standard and absolute error calculations show that for all examples Haar wavelet standard error values are mainly in the 10-1 – 10-2 range and that absolute errors depend on the transfer function in question. Accuracy of the numerical solution depend on the value of the maximum resolution level of the operational matrix, especially at sharp turns. Higher values of the parameter m provided better agreement with the analytical solution. 404 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ (a) (b) (c) (d) Fig. 4 Тransfer function 𝑋(𝑠) = 1 𝑠√𝑠 (1 − 𝑒−𝑇𝑠) inverse Laplace transform obtained analytically and numerically by Haar wavelet method (m=1024), Invlap and NILT algorithms for (a) interval [0,1) and (c) expanded interval [0,τ). Standard and absolute errors for Haar wavelet method with three different maximum resolution levels and Invlap and NILT algorithms for intervals (b) [0,1) and (d) [0,τ). 5.5. Pulse shape functions In order to investigate the application potential of numerical inversion Laplace transform using Haar wavelet operational matrix, besides examples 1-4, three common types of pulse shape functions are also presented along with a number of other inverse Laplace transforms whose analytic inverse formulas and new derived formulas are presented in Table 3. Pulse shape functions range from a square pulse that has constant value during the pulse duration, to an exponential pulse that carries a long relaxation tail, to a triangular pulse with a peak value at the half of the pulse duration. These pulses are often used in engineering and physics. Typical application is in heat transfer analysis such as transmission-line theory based photothermal modeling of composite samples where layer thermal properties are highly pulse shape function dependent [21]. Laplace transform for the square pulse is: 𝑋(𝑠) = 1 𝜏𝑠 (𝑒−𝜏1𝑠 − 𝑒−𝜏2𝑠) (33) Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 405 When we replace 1/s by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 𝑄𝐻 𝜏 ∙ (𝑒−𝜏1𝑄𝐻 −1 − 𝑒−𝜏2𝑄𝐻 −1 ) (34) then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 𝑄𝐻 𝜏 ∙ (𝑒−𝜏1𝑄𝐻 −1 − 𝑒−𝜏2𝑄𝐻 −1 ) ∙ 𝐻 (35) The analytical inverse Laplace transform of eq. (33) is: 𝑥(𝑡) = 1 𝜏 ∙ [𝐻(𝑡 − 𝜏1) − 𝐻(𝑡 − 𝜏2)] (36) Fig. 5 Comparison of different numerical inversion algorithms for square pulse, equation (36), for interval [0,1), τ=0.5, τ1=0.3 and τ2=0.8. For Haar wavelet method, maximum resolution level is m=1024. For exponential pulse: 𝑋(𝑠) = 1 𝜏2∙(𝑠+𝜏−1)2 (37) when we replace 1/s by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 𝑄𝐻 2 𝜏2∙(𝐼+𝑄𝐻∙𝜏−1)2 (38) Then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 𝑄𝐻 2 𝜏2∙(𝐼+𝑄𝐻∙𝜏−1)2 ∙ 𝐻 (39) The analytical inverse Laplace transform of eq. (37) is: 𝑥(𝑡) = 𝑡 𝜏2 ∙ 𝑒− 𝑡 𝜏 (40) 406 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ Fig. 6 Comparison of different numerical inversion algorithms for exponential pulse, equation (40), for interval [0,1) and τ=0.1. For Haar wavelet method, maximum resolution level is m=1024. Laplace transform for the triangular pulse: 𝑋(𝑠) = 2 𝑠2 ∙ (1 − 2 ∙ 𝑒−𝜏∙𝑠) (41) When we replace 1/s by Haar wavelet operational matrix QH: �̂�(𝑄𝐻) = 2𝑄𝐻 2 ∙ (𝐼 − 2 ∙ 𝑒−𝜏∙𝑄𝐻 −1 ) (42) then the inversion of Laplace transform can be calculated by: 𝑥𝑇 = [2𝑚 −2𝑚 … −2𝑚]1×𝑚 ∙ 𝐻𝑇 ∙ 2𝑄𝐻 2 ∙ (𝐼 − 2 ∙ 𝑒−𝜏∙𝑄𝐻 −1 ) ∙ 𝐻 (43) The analytical inverse Laplace transform of eq. (41) is: 𝑥(𝑡) = 2 ∙ 𝑡 − 4 ∙ 𝑡 ∙ 𝐻(𝑡 − 𝜏) + 2 ∙ 𝐻(𝑡 − 𝜏) (44) Fig. 7 Comparison of different numerical inversion algorithms for triangular pulse, equation (43), for interval [0,1) and τ=0.5. For Haar wavelet method, maximum resolution level is m=1024. Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 407 From Figs. 5-7 can be seen that numerical inversion Laplace transform using Haar wavelet operational matrix performs very well in case of all presented pulse shapes. It is in a good agreement with analytical solutions at peaks of exponential and triangular pulses. The greatest divergence from the analytical solution was seen at level changes of the square pulse where both Invlap and NILT diverged from the analytical solution as well. If we take into account all presented examples it can be seen that the Haar wavelet method and the Invlap algorithm performed better than the NILT algorithm. In order to further examine the effectiveness of Haar wavelet and Invlap numerical solutions the order of convergence was examined. For that purpose, two transfer functions were chosen, the one given in the heading 5.4 and the one for the triangular pulse given by equation 43, both of them having sharp turns at t=0.5 s. The order of convergence (OC) was calculated as [22]: 𝑂𝐶 = 𝑙𝑜𝑔 𝑥𝑖−1−𝑥𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐 𝑥𝑖−𝑥𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐 log 2 (45) where xi and xi-1 are values computed numerically in ith and (i-1)th step, respectively. In Table 1 comparison between the orders of convergence (OC) for the transfer function 𝑋(𝑠) = 1 𝑠√𝑠 (1 − 𝑒−𝑇𝑠) inverse Laplace transform obtained analytically and numerically by the Haar wavelet method and the Invlap algorithm is given. Sharp turn at t=0.5 s (Fig. 4a) was examined. From the last column of Table 1 can be seen that for the mash values of 64 and 128 the absolute error of the Haar wavelet method is greater than that of the Invlap algorithm. Higher mash values of 256 and 512 result in error ratios close to 1. Haar wavelet method has smaller absolute error than Invlap method for m=1024. However, all order of convergence values for Invlap algorithm were negative designating the reduction of accuracy. It can be concluded that, for the example in question, the Haar wavelet method performed better than Invlap, especially for higher mash values. Table 1 Comparison between the orders of convergence (OC) for the transfer function 𝑋(𝑠) = 1 𝑠√𝑠 (1 − 𝑒−𝑇𝑠) inverse Laplace transform (Fig. 4) obtained analytically and numerically by the Haar wavelet method and Invlap algorithm for t=0.5 s m Haar wavelet method Invlap algorithm Error ratio x(t) Abs. error OC x(t) Abs. error OC 64 0.7328 0.0284 0.0789 0.7021 0.0016 -5.1404 17.75 128 0.7478 0.0173 0.0547 0.7386 0.0084 -4.5483 2.06 256 0.7592 0.0097 0.0370 0.7582 0.0087 -2.0231 1.16 512 0.7680 0.0046 0.0245 0.7673 0.0039 -0.9840 1.16 1024 0.7747 0.0013 0.0160 0.7715 0.0019 -0.4885 0.68 Analytical solution: 0.7979 In Table 2 comparison between the orders of convergence (OC) for the triangular pulse inverse Laplace transform obtained analytically and numerically by the Haar wavelet method and the Invlap algorithm is given. Sharp turn at t=0.5 s (Fig. 7) was examined. From the last column of Table 2 can be seen that for the mash values of 64, 128 and 256 the absolute error of the Haar wavelet method is multiple times greater than that of the 408 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ Invlap algorithm. Mash values of 256 resulted in error ratio of 40. With higher mash values error ratio decreased. Orders of convergence for these methods were relatively low, Invlap algorithm having one order of magnitude lower OC values than Haar. For this example, Invlap algorithm performed better than Haar wavelet method. Table 2 Comparison between the orders of convergence (OC) for the triangular pulse, equation 43, inverse Laplace transform obtained analytically and numerically by Haar wavelet method and Invlap algorithm for t=0.5 s m Haar wavelet method Invlap algorithm Error ratio x(t) Abs. error OC x(t) Abs. error OC 64 0.9703 0.0140 0.4129 0.9847 0.0006 0.0144 23.33 128 0.9812 0.0109 0.3423 0.9933 0.0012 0.0732 9.08 256 0.9881 0.0080 0.2792 0.9958 0.0002 0.0745 40.00 512 0.9924 0.0056 0.2255 0.9965 0.0016 0.0468 3.50 1024 0.9952 0.0039 0.1809 0.9966 0.0024 0.0249 1.63 Analytical solution: 1.00 Table 3 Numerical inversion of Laplace transform using operational matrix QH where: 𝑥𝑇 = [2𝑚 −2𝑚 2𝑚 −2𝑚 … −2𝑚] 1×𝑚 ∙ 𝐻𝑇 ∙ �̂�(𝑄𝐻) ∙ 𝐻 X(s) x(t) �̂�(𝑸𝑯) 𝒂𝒓𝒄𝒕𝒂𝒏 𝒌 𝒔 1 𝑡 sin(𝑘 ∙ 𝑡) atan(𝑘 ∙ 𝑸𝑯) 𝟏 𝒔√𝒔 ∙ 𝒆−𝒌√𝒔 2√ 𝑡 𝜋 ∙ 𝑒− 𝑘2 4𝑡 − 𝑘 ∙ 𝑒𝑟𝑓𝑐 ( 𝑘 2√𝑡 ) 𝑸𝑯 1.5 ∙ 𝑒−𝑘𝑸𝑯 −0.5 𝒆−𝒌√𝒔 √𝒔(𝒂 + √𝒔) 𝑒𝑎𝑘 ∙ 𝑒𝑎2𝑡 ∙ 𝑒𝑟𝑓𝑐(𝑎√𝑡 + 𝑘 2√𝑡 ) 𝑸𝑯 ∙ (𝑰 + 𝑎√𝑸𝑯)−1 ∙ 𝑒−𝑘𝑄𝐻 −0.5 𝟏 𝒔√𝒔 ∙ 𝒆 𝒌 𝒔 1 √𝜋∙𝑘 sinh (2 ∙ √𝑘 ∙ 𝑡) 𝑸𝑯 1.5 ∙ 𝑒𝑘𝑸𝑯 𝟏 𝒔 ∙ 𝒆−𝒌√𝒔 𝑒𝑟𝑓𝑐 ( 𝑘 2√𝑡 ) 𝑸𝑯 ∙ 𝑒−𝑘𝑸𝑯 −0.5 𝟏 𝒔√𝒔 + 𝟏 𝑒𝑟𝑓(√𝑡) 𝑸𝑯 2 (𝑸𝑯 + 𝑸𝑯 2 ) −0.5 𝟏 √𝒔(√𝒔 + 𝒂) 𝑒𝑎2𝑡 ∙ 𝑒𝑟𝑓𝑐(𝑎√𝑡) 𝑸𝑯(𝑰 + 𝑸𝑯 0.5 ∙ 𝑎) −1 𝟏 𝒔√𝒔 ∙ 𝒆− 𝒌 𝒔 1 √𝜋 ∙ 𝑘 ∙ sin (2√𝑘𝑡) 𝑸𝑯 1.5 ∙ 𝑒−𝑘𝑸𝑯 𝟏 𝒔√𝒔 (𝟏 − 𝒆−𝑻𝒔) 2 √𝜋 ∙ √𝑥 ∙ 𝐻(𝑥) − √𝑥 − 𝑇 ∙ 𝐻(𝑥 − 𝑇) 𝑸𝑯 1.5 ∙ (𝑰 − 𝑒−𝑇𝑸𝑯 −1 ) 𝒃𝟐 − 𝒂𝟐 (𝒔 − 𝒂𝟐) ∙ (√𝒔 + 𝒃) 𝑒𝑎2 ∙ [𝑏 − 𝑎 ∙ 𝑒𝑟𝑓𝑐(𝑎√𝑡)] − −𝑏 ∙ 𝑒𝑏2𝑡 ∙ 𝑒𝑟𝑓𝑐(𝑏√𝑡) (𝑏2 − 𝑎2) ∙ 𝑸𝑯 1.5 ∙ ∙ [(𝑰 − 𝑎2 ∙ 𝑸𝑯) ∙ (𝑰 + 𝑏 ∙ 𝑸𝑯 0.5)]−1 Haar Wavelet Operational Matrix Based Numerical Inversion of Laplace Transform For Irrational ... 409 6. CONCLUSION In this study investigation of the Haar wavelet operational matrix application in the inverse Laplace transform numerical calculations for the case of irrational and transcendental transfer functions was presented. Results for a number of analytically solved inverse Laplace transforms of periodic and non-periodic functions are presented and obtained results are compared with the analytical solutions and results obtained by Invlap and NILT - algorithms that are known to be effective when irrational and transcendental functions are in question. Agreement of the numerical and analytical solutions is quantitatively evaluated using standard and absolute error calculations. For all presented examples Haar wavelet method standard error values are in the 10-1 – 10-2 range and absolute errors depend on the transfer function in question. Accuracy of the numerical solution depends on the value of the operational matrix maximum resolution level. Higher values of the operational matrix maximum resolution level improve the agreement between the numerical and analytical solutions especially at sharp turns when longer intervals require higher resolution levels. When compared to Invlap and NILT algorithms, all the algorithms used have given acceptable results. Haar wavelet method and Invlap algorithm performed better that NILT for all presented examples. In order to further examine the effectiveness of Haar wavelet and Invlap numerical solutions the order of convergence was examined in case of functions with sharp turns. Overall results showed that, although Invlap algorithm performed better than the Haar wavelet method in most cases, results obtained by Haar were in good agreement with analytical solutions for all presented examples. This approach is especially useful when the original cannot be represented by an analytical formula and numerical method must be used. In that case validity of the obtained result can be crosschecked and error can be estimated. The application potential of numerical inversion Laplace transform using Haar wavelet operational matrix is additionally confirmed on three commonly used types of pulse shape functions and a number of inverse Laplace transforms whose analytic inverse formulas and new derived formulas are presented showing that the presented method is efficient and can be easily coded. Acknowledgement: This work was financially supported by the Ministry of Science, Technological Development and Innovations of the Republic of Serbia (Contract No. 451-03-47/2023-01/200017). REFERENCES [1] G. Hariharan and K. Kannan, "A comparative study of Haar Wavelet Method and Homotropy Perturbation Method for solving one-dimensional Reaction-Diffusion Equations", Int. J. Appl. Math. Comput., vol. 3, no. 1, pp. 21-34, 2011. [2] Y. Chen, Y. Wu, Y. Cui, Z. Wang and D. Jin, "Wavelet method for a class of fractional convection-diffusion equation with variable coefficients", J. Comput. Sci., vol. 1, no. 3, pp. 146-149, 2010. [3] N. M. Bujurke, S. C. Shiralashetti and C. S. Salimath, "An application of single-term Haar wavelet series in the solution of nonlinear oscillator equations", J. Comput. Appl. Math., vol. 227, no. 2, pp. 234-244, 2009. [4] V. Bruni, B. Piccoli and D. Vitulano, "Wavelets and partial differential equations for image denoising", Electron. Lett. Comput. Vis. Image Anal., vol. 6, no. 2, pp. 36-53, 2008. [5] G. Hariharan and K. Kannan, "An Overview of Haar Wavelet Method for Solving Differential and Integral Equations", World Appl. Sci. J., vol. 23, no. 12, pp. 1-14, 2013. [6] C. H. Hsiao and S.P. Wu, "Numerical solution of time-varying functional differential equations via Haar wavelets", Appl. Math. Comput., vol. 188, no. 1, pp. 1049-1058, 2007. [7] N. Salamat, M. M. S. Missen and V. B. Surya Prasath, "Recent developments in computational color image denoising with PDEs to deep learning: a review", Artif. Intell. Rev., vol. 54, pp. 6245–6276, 2021. 410 Z. STANIMIROVIĆ, I. STANIMIROVIĆ, S. GALOVIĆ, K. ĐORĐEVIĆ, E. SULJOVRUJIĆ [8] C. Tian and Y. Chen, "Image Segmentation and Denoising Algorithm Based on Partial Differential Equations", IEEE Sensors J., vol. 20, no. 20, pp. 11935-11942, 2020. [9] W. Cai, W. Chen, J. Fang and S. Holm, "A Survey on Fractional Derivative Modeling of Power-Law Frequency- Dependent Viscous Dissipative and Scattering Attenuation in Acoustic Wave Propagation", ASME. Appl. Mech. Rev., vol. 70, no. 3, p. 030802, 2018. [10] T. Rabczuk, H. Ren and X. Zhuang, "A Nonlocal Operator Method for Partial Differential Equations with Application to Electromagnetic Waveguide Problem", Comput. Mater. Continua, vol. 59, no. 1, pp. 31-55, 2019. [11] H. Sheng, Y. Li and Y.Q. Chen, "Application of numerical inverse Laplace transform algorithms in fractional calculus", J. Franklin Institute, vol. 348, no. 2, pp. 315-330, 2011. [12] S. M. Aznam, A. Hussin, "Numerical Method for Inverse Laplace Transform with Haar Wavelet Operational Matrix", Mal. J. Fundam. Appl. Sci., vol. 8, no. 4, pp. 204-210, 2012. [13] J. L. Wu, C. H. Chen and C. F. Chen, "Numerical inversion of Laplace transform using Haar wavelet operational matrices", IEEE Trans. Circ. Syst. I: Fundam. Theory Appl., vol. 48, no. 1, pp. 120-122, 2001. [14] C. H. Hsiao, "Numerical inversion of Laplace transform via wavelet in ordinary differential equations", Comput. Methods Differ. Equ., vol. 3, pp. 186-194, 2014. [15] J. A. C. Weideman, "Algorithms for Parameter Selection in the Weeks Method for Inverting the Laplace Transform", SIAM J. Sci. Comput., vol. 21, pp. 111-128, 1999. [16] J. Kotyk, "MATLAB Algorithms for The Laplace Transform Inversion", In Proceedings of the MATLAB Conference, Prague, Czech Republic, 2005, pp. 1-19. [17] F. R. De Hoog, J. H. Knight and A. N. Stokes, "An improved method for numerical inversion of Laplace transforms", SIAM J. Sci. Stat. Comput., vol. 3, no. 3, pp. 357-366, 1982. [18] H. Stehfest, "Algorithm 368: Numerical Inversion of Laplace Transform", Commun. ACM, vol. 13, no. 1, pp. 47- 49, 1970. [19] N. Al-Zubaidi Smith and L. Brancik, "Comparative Study on One-Dimensional Numerical Inverse Laplace Transform Methods for Electrical Engineering", Elektrorevue, vol. 18, no. 1, pp. 1-7, 2016. [20] D. J. Halsted and D. E. Brown, "Zakian's technique for inverting Laplace transforms", Chem. Eng. J., vol. 3, pp. 312-313, 1972. [21] G. Chen and P. Huib, "Pulsed photothermal modeling of composite samples based on transmission-line theory of heat conduction", Thin Solid Films, vol. 339, pp. 58-67, 1999. [22] M. Sorrenti, M. Di Sciuva, J. Majak and F. Auriemma, "Static Response and Buckling Loads of Multilayered Composite Beams Using the Refined Zigzag Theory and Higher-Order Haar Wavelet Method", Mech. Compos. Mater., vol. 57, no.1, 2021.