An Efficient Exponential Sum Approximation of Power-Law Kernels for Solving Fractional Differential Equation

TL;DR

This paper introduces a computationally efficient exponential sum approximation for power-law kernels used in fractional differential equations, enabling faster and accurate solutions with reduced computational cost.

Contribution

It proposes a novel finite exponential sum approximation method with error minimization and a two-stage compression algorithm based on Prony's method for solving fractional differential equations.

Findings

The method achieves high accuracy in approximating power-law kernels.

Numerical results demonstrate the efficiency and convergence of the proposed approach.

The approach reduces computational complexity compared to existing methods.

Abstract

In this work, we present a comprehensive framework for approximating the weakly singular power-law kernel $t^{\alpha-1}$ of fractional integral and differential operators, where $\alpha \in (0,1)$ and $t \in [\delta,T]$ with $0<\delta<T<\infty$, using a finite sum of exponentials. This approximation method begins by substituting an exponential function into the Laplace transform of the power function, followed by the application of the trapezoidal rule to approximate the resulting integral. To ensure computational feasibility, the integral limits are truncated, leading to a finite exponential sum representation of the kernel. In contrast to earlier approaches, we pre-specify the admitted computational cost (measured in terms of the number of exponentials) and minimize the approximation error. Furthermore, to reduce the computational cost while maintaining accuracy, we present a two-stage algorithm based on Prony's method that compresses the exponential sum. The compressed kernel is then embedded into the Riemann-Liouville fractional integral and applied to solve fractional differential equations. To this end, we discuss two solution strategies, namely (a) method based on piecewise constant interpolation and (b) a transformation of the original fractional differential equation into a system of first-order ordinary differential equations (ODEs). This reformulation makes the problem solvable by standard ODE solvers with low computational cost while retaining the accuracy benefits of the exponential-sum-approximation. Finally, we apply the proposed strategies to solve some well-known fractional differential equations and demonstrate the advantages, accuracy, and the experimental order of convergence of the methods through numerical results.

Tables

Table 1. Table 1: Max. error values obtained before and after the application of Prony’s method for various values of L L and α ∈ { 0.1 , 0.5 , 0.9 } \alpha\in\{0.1,0.5,0.9\} over the interval [ 10 − 2 , 1 ] [10^{-2},1] . In all cases, the numbers of non-positive nodes M M , the values of L p L_{\text{p}} and K K in Prony’s method and the final number of exponential terms after reduction L f L_{\text{f}} are provided. Here, the value of ϵ \epsilon is set to 10 − 10 10^{-10} .

$𝜶$	$𝑳$	$𝑴$	${𝑳}_{\text{p}}$	$𝑲$	${𝑳}_{\text{f}}$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},\mathrm{𝟏}]}\left|𝒆(𝒕,{𝜽}_{\mathrm{𝟎}})\right|$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},\mathrm{𝟏}]}\left|{𝒆}_{\text{p}}\left(𝒕\right)\right|$
$0.1$	$32$	24	24	1	9	5.102467$\left(-2\right)$	5.102386$\left(-2\right)$
	$64$	49	49	3	18	6.483821$\left(-6\right)$	6.510213$\left(-6\right)$
	$128$	98	98	4	34	1.320726$\left(-8\right)$	1.980379$\left(-10\right)$
	$256$	196	196	4	64	7.207191$\left(-9\right)$	6.868319$\left(-10\right)$
$0.5$	$32$	27	27	1	6	8.401490$\left(-2\right)$	8.401582$\left(-2\right)$
	$64$	55	55	2	11	3.577202$\left(-4\right)$	3.577193$\left(-4\right)$
	$128$	110	110	4	22	3.988015$\left(-9\right)$	3.802676$\left(-9\right)$
	$256$	220	220	5	41	3.518998$\left(-10\right)$	5.593037$\left(-11\right)$
$0.9$	$128$	124	124	2	6	4.027975$\left(-3\right)$	4.027975$\left(-3\right)$
	$256$	248	248	2	10	2.591330$\left(-5\right)$	2.591330$\left(-5\right)$
	$512$	496	496	4	20	1.240738$\left(-9\right)$	1.240076$\left(-9\right)$
	$1024$	993	993	5	36	1.342770$\left(-11\right)$	1.039657$\left(-11\right)$

Table 2. Table 2: Errors and EOC at T = 1 T=1 for ( 49 ) with α = 0.1 \alpha=0.1 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	8.08(-2)		7.89(-2)		2.30(-2)		2.05(-2)		2.14(-4)		1.73(-3)
$2^{-3}$	4.77(-2)	0.76	4.57(-2)	0.79	1.78(-2)	0.37	1.55(-2)	0.40	1.71(-4)	0.32	2.13(-3)	***
$2^{-4}$	2.80(-2)	0.77	2.60(-2)	0.81	1.17(-2)	0.60	9.60(-3)	0.70	1.10(-3)	***	8.61(-4)	1.31
$2^{-5}$	1.65(-2)	0.76	1.45(-2)	0.84	7.21(-3)	0.70	5.14(-3)	0.90	1.70(-3)	***	2.63(-4)	1.71
$2^{-6}$	9.95(-3)	0.73	7.98(-3)	0.87	4.60(-3)	0.65	2.58(-3)	1.00	1.89(-3)	***	7.00(-5)	1.91
$2^{-7}$	6.29(-3)	0.66	4.31(-3)	0.89	3.26(-3)	0.50	1.26(-3)	1.03	1.95(-3)	***	1.68(-5)	2.06
$2^{-8}$	4.28(-3)	0.56	2.31(-3)	0.90	2.60(-3)	0.33	6.17(-4)	1.03	1.96(-3)	***	3.55(-6)	2.24
$2^{-9}$	3.19(-3)	0.42	1.22(-3)	0.92	2.27(-3)	0.19	3.02(-4)	1.03	1.96(-3)	***	5.67(-7)	2.65
$2^{-10}$	2.61(-3)	0.29	6.45(-4)	0.93	2.12(-3)	0.10	1.48(-4)	1.02	1.96(-3)	***	4.58(-9)	6.95
${𝑳}_{\text{f}}$	24		50		24		50		24		50

Table 3. Table 3: Errors and EOC at T = 1 T=1 for ( 49 ) with α = 0.5 \alpha=0.5 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	3.30(-1)		3.30(-1)		1.01(-1)		1.01(-1)		3.34(-2)		3.34(-2)
$2^{-3}$	1.51(-1)	1.12	1.51(-1)	1.12	8.48(-2)	0.26	8.47(-2)	0.26	1.38(-2)	1.27	1.38(-2)	1.27
$2^{-4}$	7.55(-2)	1.01	7.55(-2)	1.01	4.86(-2)	0.80	4.85(-2)	0.80	3.96(-3)	1.80	3.98(-3)	1.80
$2^{-5}$	3.78(-2)	1.00	3.78(-2)	1.00	2.42(-2)	1.01	2.41(-2)	1.01	9.73(-4)	2.03	9.95(-4)	2.00
$2^{-6}$	1.89(-2)	1.00	1.89(-2)	1.00	1.16(-2)	1.06	1.15(-2)	1.06	2.08(-4)	2.22	2.32(-4)	2.10
$2^{-7}$	9.50(-3)	1.00	9.47(-3)	1.00	5.56(-3)	1.06	5.54(-3)	1.06	3.10(-5)	2.75	5.22(-5)	2.15
$2^{-8}$	4.75(-3)	1.00	4.73(-3)	1.00	2.71(-3)	1.04	2.68(-3)	1.04	9.71(-6)	1.68	1.14(-5)	2.19
$2^{-9}$	2.38(-3)	0.99	2.36(-3)	1.00	1.34(-3)	1.01	1.32(-3)	1.03	1.99(-5)	***	2.44(-6)	2.23
$2^{-10}$	1.20(-3)	0.99	1.18(-3)	1.00	6.74(-4)	0.99	6.52(-4)	1.02	2.13(-5)	***	4.78(-7)	2.35
${𝑳}_{\text{f}}$	17		33		17		33		17		33

Table 4. Table 4: Errors and EOC at T = 1 T=1 for ( 49 ) with α = 0.9 \alpha=0.9 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟐𝟓𝟔}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟐𝟓𝟔}$		$𝑳=\mathrm{𝟐𝟓𝟔}$		$𝑳=\mathrm{𝟓𝟏𝟐}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	4.22(-1)		4.25(-1)		2.64(-1)		2.66(-1)		6.11(-2)		6.11(-2)
$2^{-3}$	1.98(-1)	1.09	1.99(-1)	1.09	1.94(-1)	0.44	1.95(-1)	0.44	1.83(-2)	1.74	1.83(-2)	1.74
$2^{-4}$	9.86(-2)	1.01	9.99(-2)	1.00	1.03(-1)	0.91	1.04(-1)	0.91	4.71(-3)	1.96	4.70(-3)	1.96
$2^{-5}$	4.86(-2)	1.02	5.02(-2)	0.99	4.94(-2)	1.06	5.08(-2)	1.04	1.16(-3)	2.01	1.16(-3)	2.02
$2^{-6}$	2.36(-2)	1.04	2.52(-2)	0.99	2.30(-2)	1.10	2.45(-2)	1.05	2.89(-4)	2.01	2.83(-4)	2.03
$2^{-7}$	1.11(-2)	1.09	1.26(-2)	1.00	1.05(-2)	1.13	1.19(-2)	1.04	7.53(-5)	1.94	6.94(-5)	2.03
$2^{-8}$	4.80(-3)	1.22	6.31(-3)	1.00	4.39(-3)	1.26	5.89(-3)	1.02	2.31(-5)	1.70	1.70(-5)	2.02
$2^{-9}$	1.64(-3)	1.55	3.15(-3)	1.00	1.41(-3)	1.63	2.91(-3)	1.01	1.02(-5)	1.17	4.22(-6)	2.02
$2^{-10}$	7.42(-5)	4.47	1.57(-3)	1.00	4.67(-5)	4.92	1.44(-3)	1.01	7.06(-6)	0.54	1.04(-6)	2.01
${𝑳}_{\text{f}}$	9		17		9		17		17		32

Table 5. Table 5: Errors and EOC at T = 10 T=10 for ( 50 ) with α = 0.1 \alpha=0.1 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟐𝟓𝟔}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	2.58(-4)		1.50(-4)		4.74(-4)		3.59(-4)		2.35(-4)		2.35(-4)
$2^{-3}$	1.90(-4)	0.44	8.20(-5)	0.87	2.91(-4)	0.70	1.79(-4)	1.01	1.08(-4)	1.12	1.08(-4)	1.12
$2^{-4}$	1.53(-4)	0.32	4.41(-5)	0.89	1.99(-4)	0.54	8.89(-5)	1.01	5.04(-5)	1.11	5.04(-5)	1.11
$2^{-5}$	1.32(-4)	0.21	2.35(-5)	0.91	1.54(-4)	0.37	4.41(-5)	1.01	2.34(-5)	1.10	2.34(-5)	1.10
$2^{-6}$	1.21(-4)	0.12	1.24(-5)	0.92	1.31(-4)	0.23	2.19(-5)	1.01	1.09(-5)	1.10	1.09(-5)	1.10
$2^{-7}$	1.15(-4)	0.07	6.55(-6)	0.93	1.20(-4)	0.13	1.08(-5)	1.01	5.08(-6)	1.10	5.08(-6)	1.10
$2^{-8}$	1.12(-4)	0.04	3.42(-6)	0.94	1.14(-4)	0.07	5.39(-6)	1.01	2.37(-6)	1.10	2.37(-6)	1.10
$2^{-9}$	1.11(-4)	0.02	1.78(-6)	0.94	1.12(-4)	0.04	2.68(-6)	1.01	1.10(-6)	1.10	1.10(-6)	1.10
$2^{-10}$	1.10(-4)	0.01	9.24(-7)	0.95	1.10(-4)	0.02	1.33(-6)	1.01	5.15(-7)	1.10	5.15(-7)	1.10
${𝑳}_{\text{f}}$	27		54		27		54		54		105

Table 6. Table 6: Errors and EOC at T = 10 T=10 for ( 50 ) with α = 0.5 \alpha=0.5 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟔𝟒}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟏𝟐𝟖}$		$𝑳=\mathrm{𝟐𝟓𝟔}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	4.14(-3)		1.69(-3)		4.33(-3)		1.77(-3)		2.01(-4)		2.01(-4)
$2^{-3}$	3.34(-3)	0.31	8.52(-4)	0.99	3.37(-3)	0.36	8.19(-4)	1.11	6.91(-5)	1.54	6.91(-5)	1.54
$2^{-4}$	2.94(-3)	0.18	4.28(-4)	0.99	2.93(-3)	0.20	3.87(-4)	1.08	2.39(-5)	1.53	2.39(-5)	1.53
$2^{-5}$	2.74(-3)	0.10	2.14(-4)	0.99	2.73(-3)	0.10	1.86(-4)	1.06	8.36(-6)	1.52	8.36(-6)	1.52
$2^{-6}$	2.64(-3)	0.05	1.07(-4)	1.00	2.63(-3)	0.05	9.07(-5)	1.04	2.92(-6)	1.51	2.93(-6)	1.51
$2^{-7}$	2.58(-3)	0.03	5.39(-5)	1.00	2.58(-3)	0.03	4.45(-5)	1.03	1.02(-6)	1.51	1.03(-6)	1.51
$2^{-8}$	2.56(-3)	0.01	2.70(-5)	1.00	2.56(-3)	0.01	2.19(-5)	1.02	3.58(-7)	1.52	3.63(-7)	1.51
$2^{-9}$	2.55(-3)	0.01	1.35(-5)	1.00	2.55(-3)	0.01	1.08(-5)	1.01	1.23(-7)	1.54	1.28(-7)	1.50
$2^{-10}$	2.54(-3)	0.01	6.77(-6)	1.00	2.54(-3)	0.01	5.39(-6)	1.01	4.05(-8)	1.61	4.51(-8)	1.50
${𝑳}_{\text{f}}$	18		37		18		37		37		72

Table 7. Table 7: Errors and EOC at T = 10 T=10 for ( 50 ) with α = 0.9 \alpha=0.9 obtained by different methods.

	CI				BE				TR
	$𝑳=\mathrm{𝟐𝟓𝟔}$		$𝑳=\mathrm{𝟓𝟏𝟐}$		$𝑳=\mathrm{𝟐𝟓𝟔}$		$𝑳=\mathrm{𝟓𝟏𝟐}$		$𝑳=\mathrm{𝟓𝟏𝟐}$		$𝑳=\mathrm{𝟏𝟎𝟐𝟒}$
$𝒉$	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC	Error	EOC
$2^{-2}$	7.75(-4)		7.70(-4)		8.47(-4)		8.44(-4)		1.47(-5)		1.47(-5)0
$2^{-3}$	3.76(-4)	1.04	3.71(-4)	1.05	3.80(-4)	1.15	3.77(-4)	1.16	3.78(-6)	1.97	3.78(-6)0	1.97
$2^{-4}$	1.87(-4)	1.01	1.82(-4)	1.03	1.82(-4)	1.06	1.78(-4)	1.08	9.68(-7)	1.97	9.68(-7)0	1.97
$2^{-5}$	9.56(-5)	0.97	9.04(-5)	1.01	9.14(-5)	1.00	8.65(-5)	1.04	2.47(-7)	1.97	2.47(-7)0	1.97
$2^{-6}$	5.02(-5)	0.93	4.50(-5)	1.01	4.77(-5)	0.94	4.26(-5)	1.02	6.32(-8)	1.97	6.34(-8)0	1.97
$2^{-7}$	2.77(-5)	0.86	2.24(-5)	1.00	2.63(-5)	0.86	2.11(-5)	1.01	1.60(-8)	1.98	1.62(-8)0	1.97
$2^{-8}$	1.65(-5)	0.75	1.12(-5)	1.00	1.57(-5)	0.74	1.05(-5)	1.01	3.98(-9)	2.01	4.16(-9)0	1.96
$2^{-9}$	1.09(-5)	0.60	5.60(-6)	1.00	1.05(-5)	0.58	5.25(-6)	1.00	8.85(-10)	2.17	1.06(-9)0	1.96
$2^{-10}$	8.09(-6)	0.43	2.80(-6)	1.00	7.91(-6)	0.41	2.62(-6)	1.00	9.16(-11)	3.27	2.74(-10)	1.96
${𝑳}_{\text{f}}$	19		37		19		37		37		70

Table 8. Table 8: Max. error values obtained before and after the application of Prony’s method for various values of L L and α ∈ { 0.1 , 0.5 , 0.9 } \alpha\in\{0.1,0.5,0.9\} over the interval [ 10 − 2 , 10 3 ] [10^{-2},10^{3}] . In all cases, the numbers of non-positive nodes M M , the values of L p L_{\text{p}} and K K in Prony’s method and the final number of exponential terms after reduction L f L_{\text{f}} are provided. Here, the value of ϵ \epsilon is set to 10 − 10 10^{-10} .

$𝜶$	$𝑳$	$𝑴$	${𝑳}_{\text{p}}$	$𝑲$	${𝑳}_{\text{f}}$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}]}\left|𝒆(𝒕,{𝜽}_{\mathrm{𝟎}})\right|$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}]}\left|{𝒆}_{\text{p}}\left(𝒕\right)\right|$
$0.1$	$32$	20	20	1	13	2.472386(-1)0	2.472386(-1)0
	$64$	41	41	1	24	1.335328(-4)0	1.335343(-4)0
	$128$	81	81	3	50	1.582185(-8)0	1.542839(-10)
	$256$	163	163	3	96	8.214528(-9)0	6.045368(-10)
$0.5$	$32$	24	24	1	9	2.050205(-1)0	2.050205(-1)0
	$64$	49	49	2	17	1.159256(-3)0	1.159256(-3)0
	$128$	98	98	3	33	4.535111(-8)0	4.532144(-8)0
	$256$	195	195	4	65	3.326726(-10)	4.989634(-12)
$0.9$	$128$	121	121	2	9	4.834038(-3)0	4.834038(-3)0
	$256$	242	242	3	17	3.344168(-5)0	3.344168(-5)0
	$512$	484	484	4	32	2.046576(-9)0	2.046172(-9)0
	$1024$	968	968	5	61	8.260503(-12)	5.237699(-12)

Table 9. Table 9: Max. error values obtained before and after the application of Prony’s method for fixed values of L L and α ∈ { 0.1 , 0.5 , 0.9 } \alpha\in\{0.1,0.5,0.9\} over the interval [ 10 − 2 , 10 3 ] [10^{-2},10^{3}] . In all cases, the numbers of non-popsitive nodes, M M , the values of L p L_{\text{p}} and K K in Prony’s method, and the final number of exponential terms after reduction L f L_{\text{f}} are provided.

$𝜶$	$ϵ$	$𝑴$	${𝑳}_{\text{p}}$	$𝑲$	${𝑳}_{\text{f}}$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}]}\left|𝒆(𝒕,{𝜽}_{\mathrm{𝟎}})\right|$	${\mathrm{𝐦𝐚𝐱}}_{𝒕\in [{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}]}\left|{𝒆}_{\text{p}}\left(𝒕\right)\right|$
$0.1$ $\left(L=256\right)$	${10}^{-10}$	163	163	3	96	8.214528(-9)0	6.045368(-10)
	${10}^{-11}$	168	168	3	91	9.331558(-10)	5.064520(-10)
	${10}^{-12}$	173	173	4	87	1.058709(-10)	1.094236(-12)
	${10}^{-13}$	177	177	4	83	1.196554(-11)	6.685971(-13)
$0.5$ $\left(L=256\right)$	${10}^{-10}$	195	195	4	65	3.326726(-10)	4.989634(-12)
	${10}^{-11}$	199	199	4	61	3.741008(-11)	3.887078(-12)
	${10}^{-12}$	203	203	5	58	4.185097(-12)	3.197442(-14)
	${10}^{-13}$	206	206	5	55	5.879741(-13)	1.278977(-13)
$0.9$ $\left(L=1024\right)$	${10}^{-10}$	968	968	5	61	8.260503(-12)	5.237699(-12)
	${10}^{-11}$	972	972	5	57	8.537615(-13)	5.442313(-13)
	${10}^{-12}$	975	975	5	54	8.926193(-14)	5.584422(-14)
	${10}^{-13}$	978	978	5	51	3.264056(-14)	2.842171(-14)

Table 10. Table 10: Values of the weights ρ k \rho_{k} and exponents η k \eta_{k} obtained by Prony’s method for various values of α ∈ { 0.1 , 0.5 , 0.9 } \alpha\in\{0.1,0.5,0.9\} over the intervals [ 10 − 2 , 1 ] [10^{-2},1] , [ 10 − 5 , 1 ] [10^{-5},1] , and [ 10 − 2 , 10 3 ] [10^{-2},10^{3}] . Here, the value of ϵ \epsilon is set to 10 − 10 10^{-10} .

	$[{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},\mathrm{𝟏}]$		$[{\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}},\mathrm{𝟏}]$		$[{\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}},{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}]$
$𝜶$	${𝝆}_{𝒌}$	${𝜼}_{𝒌}$	${𝝆}_{𝒌}$	${𝜼}_{𝒌}$	${𝝆}_{𝒌}$	${𝜼}_{𝒌}$
$0.1$ $\left(L=256\right)$	0.1887	$-0.8580$	0.3076	$-0.8810$	$0.6138\left(-3\right)$	$-0.8810\left(-3\right)$
	0.3202	$-0.6074$	0.4804	$-0.4834$	$0.9586\left(-3\right)$	$-0.4834\left(-3\right)$
	0.3384	$-0.2926$	0.3347	$-0.1006$	$0.6677\left(-3\right)$	$-0.1006\left(-3\right)$
	0.2033	$-0.0569$
$0.5$ $\left(L=256\right)$	0.2239	$-0.9500$	0.2600	$-0.8070$	0.0082	$-0.8070\left(-3\right)$
	0.3026	$-0.7184$	0.4149	$-0.5384$	0.0131	$-0.5384\left(-3\right)$
	0.4290	$-0.4413$	0.5775	$-0.2336$	0.0183	$-0.2336\left(-3\right)$
	0.5265	$-0.1795$	0.6670	$-0.0284$	0.0211	$-0.0284\left(-3\right)$
	0.5778	$-0.0212$
$0.9$ $\left(L=1024\right)$	0.2618	$-0.8433$	0.2715	$-0.9807$	0.1361	$-0.9807\left(-3\right)$
	0.3766	$-0.6154$	0.3836	$-0.7131$	01922	$-0.7131\left(-3\right)$
	0.6551	$-0.3579$	0.6648	$-0.4145$	0.3332	$-0.4145\left(-3\right)$
	1.2235	$-0.1258$	1.2416	$-0.1458$	0.6223	$-0.1458\left(-3\right)$
	7.4423	$-0.0034$	7.5525	$-0.0040$	3.7852	$-0.0040\left(-3\right)$

Equations

$$w(t)\approx \sum _{l=1}^P{b}_l{\mathrm{e}}^{-a_{l}t},$$

$$t^{\alpha -1}=\mathcal{L}\{{\xi }^{-\alpha }/\Gamma (1-\alpha )\}(t),$$

$$I^{\alpha }f(t)=\frac{1}{\Gamma (\alpha )}{\int }_0^t(t-\tau {)}^{\alpha -1}f(\tau )\mathrm{d}\tau .$$

$$t^{\alpha -1}=\frac{1}{\Gamma (1-\alpha )}{\int }_0^{\infty }{\mathrm{e}}^{-pt}{p}^{-\alpha }\mathrm{d}p.$$

$$t^{\alpha -1}=\frac{1}{\Gamma (1-\alpha )}{\int }_{-\infty }^{\infty }{\mathrm{e}}^{(1-\alpha )\omega }{\mathrm{e}}^{-t{\mathrm{e}}^{\omega }}\mathrm{d}\omega .$$

$$t^{\alpha-1}\approx\frac{1}{\Gamma(1-\alpha)}\int_{l_{\min}}^{l_{\max}}g(\omega)\,\mathrm{d}\omega\approx\frac{h}{\Gamma(1-\alpha)}\bigg{(}\frac{g(\omega_{1})+g(\omega_{L})}{2}+\sum_{l=2}^{L-1}g(\omega_{l})\bigg{)}$$

$$g(\omega )={\mathrm{e}}^{(1-\alpha )\omega }{\mathrm{e}}^{-t{\mathrm{e}}^{\omega }}.$$

$$t^{\alpha -1}\approx f(t,\theta ):=\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^L{w}_l{\mathrm{e}}^{b_{l}t}\text{for}t\in [\delta ,T],$$

$$w_1$$

$$w_l$$

$$b_l$$

$$e(t,\theta )=t^{\alpha -1}-\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^L{w}_l{\mathrm{e}}^{b_{l}t}$$

$${\tilde{w}}_l=T^{1-\alpha }{w}_l,{\tilde{b}}_l=T{b}_l,(l=1,\dots ,L)$$

$$e(t,\theta )=T^{\alpha -1}e(\tilde{t},\tilde{\theta }).$$

$$I_{\min }={\int }_{-\infty }^{l_{\min }}g(\omega )\mathrm{d}\omega \text{and}{I}_{\max }={\int }_{l_{\max }}^{\infty }g(\omega )\mathrm{d}\omega$$

$$I_{\min }=\frac{1}{t^{1-\alpha }}{\int }_0^{t{\mathrm{e}}^{l_{\min }}}{v}^{-\alpha }{\mathrm{e}}^{-v}\mathrm{d}v=\frac{1}{t^{1-\alpha }}\gamma (1-\alpha ,t{\mathrm{e}}^{l_{\min }})$$

$$\gamma (\beta ,q)={\int }_0^q{p}^{\beta -1}{\mathrm{e}}^{-p}\mathrm{d}p\text{for}\beta >0\text{and}q>0,$$

$$l_{\min }\le \frac{\ln (ϵ(1-\alpha ))}{1-\alpha }.$$

$$l_{\min }\le \ln \left(\frac{ϵ}{T}\right).$$

$$l_{\min }=\min \left\{\ln \left(\frac{ϵ}{T}\right),\frac{\ln (ϵ(1-\alpha ))}{1-\alpha }\right\}.$$

$$I_{\max }=\frac{1}{t^{1-\alpha }}{\int }_{t{\mathrm{e}}^{l_{\max }}}^{\infty }{v}^{-\alpha }{\mathrm{e}}^{-v}\mathrm{d}v=\frac{1}{t^{1-\alpha }}\Gamma (1-\alpha ,t{\mathrm{e}}^{l_{\max }}),$$

$$\Gamma (\beta ,q)={\int }_q^{\infty }{p}^{\beta -1}{\mathrm{e}}^{-p}\mathrm{d}p\text{for}\beta >0\text{ and }q>0.$$

$$I_{\max }\approx \frac{1}{t^{1-\alpha }}(t{\mathrm{e}}^{l_{\max }}{)}^{-\alpha }\exp (-t{\mathrm{e}}^{l_{\max }})=\frac{1}{t}{\mathrm{e}}^{-\alpha l_{\max }}\exp (-t{\mathrm{e}}^{l_{\max }}).$$

$$l_{\max }(t)\approx \ln \left(-\frac{\ln ϵ}{t}\right)\text{for}\delta \le t\le T.$$

$$l_{\max }=\ln \left(-\frac{\ln ϵ}{\delta }\right).$$

$$\sum _{l=1}^L{w}_l{\mathrm{e}}^{b_{l}t}=\sum _{{\omega }_l\le 0}{w}_l{\mathrm{e}}^{b_{l}t}+\sum _{{\omega }_l>0}{w}_l{\mathrm{e}}^{b_{l}t},$$

$$\underset{l=1,2,\dots ,256}{\max }{w}_l=\{\begin{array}{ll}123.7369& \text{for }\alpha =0.1,\\ 9.3188& \text{for }\alpha =0.5,\text{ and}\\ 2.0041& \text{for }\alpha =\mathrm{0.9.}\end{array}$$

$$M:=\max \{l\in 1,2,\dots ,L:{\omega }_l\le 0\}$$

$$\sum _{k=1}^K{\rho }_k{\mathrm{e}}^{{\eta }_{k}t},$$

$${\epsilon }_{\text{p}}(t)=\frac{1}{\Gamma (1-\alpha )}\left(\sum _{l=1}^{L_{\text{p}}}{w}_l{\mathrm{e}}^{b_{l}t}-\sum _{k=1}^K{\rho }_k{\mathrm{e}}^{{\eta }_{k}t}\right)$$

$$L_{\text{f}}=K+L-L_{\text{p}}.$$

Full text

An Efficient Exponential Sum Approximation of Power-Law Kernels for Solving Fractional Differential Equations111This work was supported by the German Federal Ministry of Research, Technology and Space under Grant No. 05M22WHA.

Renu Chaudhary

Faculty of Applied Natural Sciences and Humanities (FANG), Technical University of Applied Sciences Würzburg-Schweinfurt, Ignaz-Schön-Str. 11, 97421 Schweinfurt, Germany

Kai Diethelm

Faculty of Applied Natural Sciences and Humanities (FANG), Technical University of Applied Sciences Würzburg-Schweinfurt, Ignaz-Schön-Str. 11, 97421 Schweinfurt, Germany

Afshin Farhadi Corresponding author

Email addresses: renu.chaudhary@thws.de (R. Chaudhary), kai.diethelm@thws.de (K. Diethelm),

afshin.farhadi@thws.de (A. Farhadi), fredalbert.fuchs@thws.de (F.A. Fuchs). Faculty of Applied Natural Sciences and Humanities (FANG), Technical University of Applied Sciences Würzburg-Schweinfurt, Ignaz-Schön-Str. 11, 97421 Schweinfurt, Germany

Fred A. Fuchs

Faculty of Applied Natural Sciences and Humanities (FANG), Technical University of Applied Sciences Würzburg-Schweinfurt, Ignaz-Schön-Str. 11, 97421 Schweinfurt, Germany

 

Abstract. In this work, we present a comprehensive framework for approximating the weakly singular power-law kernel $t^{\alpha-1}$ of fractional integral and differential operators, where $\alpha\in(0,1)$ and $t\in[\delta,T]$ with $0<\delta<T<\infty$, using a finite sum of exponentials. This approximation method begins by substituting an exponential function into the Laplace transform of the power function, followed by the application of the trapezoidal rule to approximate the resulting integral. To ensure computational feasibility, the integral limits are truncated, leading to a finite exponential sum representation of the kernel. In contrast to earlier approaches, we pre-specify the admitted computational cost (measured in terms of the number of exponentials) and minimize the approximation error. Furthermore, to reduce computational cost while maintaining accuracy, we present a two-stage algorithm based on Prony’s method that compresses the exponential sum. The compressed kernel is then embedded into the Riemann-Liouville fractional integral and applied to solve fractional differential equations. To this end, we discuss two solution strategies, namely (a) method based on piecewise constant interpolation and (b) a transformation of the original fractional differential equation into a system of first-order ordinary differential equations (ODEs). This reformulation makes the problem solvable by standard ODE solvers with low computational cost while retaining the accuracy benefits of the exponential-sum approximation. Finally, we apply the proposed strategies to solve some well-known fractional differential equations and demonstrate the advantages, accuracy and the experimental order of convergence of the methods through numerical results.

Keywords. Fractional differential equations, power-law kernel, exponential sum approximation, trapezoidal rule, Prony’s method

1. Introduction

The Riemann-Liouville fractional integral is a fundamental concept in fractional calculus, characterized (up to a multiplicative constant) by having the form of a Laplace-type convolution with kernel function $w(t)=t^{\alpha-1}$ where $\alpha>0$ and $t\in(0,T]$. Therefore, the accurate and efficient approximation of power-law kernels plays a crucial role in the computation of fractional integrals as required, e.g., in typical algorithms for solving fractional differential equations. In the vast majority of applications [6, 7, 8, 10, 16, 21, 22], one needs to deal with the case that $0<\alpha<1$ . In this case, the singularity at $t=0$ and the slow decay of the kernel make direct numerical computations challenging. The sum-of-exponentials method addresses this challenge by approximating the kernel $t^{\alpha-1}$ with a sum of exponential functions which are easier to handle numerically and lead to more efficient algorithms.

A variety of sum-of-exponentials methods has been introduced for this purpose. For example, Baffet [2] developed a composite Gauss–Jacobi quadrature on the Laplace-domain integral representation of $t^{\alpha-1}$. By subdividing the integral into finite intervals and truncating the infinite tails, he obtained a rapidly convergent approximation of the form

[TABLE]

valid for $t\in[\delta,T]$, with real coefficients $a_{l},b_{l}$ chosen by the quadrature. This scheme delivers high accuracy with a modest number of terms $P$, uniformly in $\alpha\in(0,1)$, and provides explicit error bounds. He also incorporated the compressed kernel into a time-stepping algorithm to solve fractional differential equations. Baffet and Hesthaven [3] split the fractional integral of a function $f$ into local and history parts and proposed a kernel compression scheme based on a multipole approximation of the Laplace transform of the kernel function. This method transforms the history part into a linear combination of auxiliary variables, each defined by a standard ordinary differential equation, so that the fractional problem can be converted to a system of ODEs. They also derived estimates for the number of auxiliary variables required to attain a prescribed error tolerance. In a different approach, Li [13] exploited the representation of the kernel as a Laplace transform,

[TABLE]

and constructed a fixed $Q$-point quadrature rule for the integral $\int_{0}^{\infty}\mathrm{e}^{-\xi t}\xi^{-\alpha}\,\mathrm{d}\xi$ that is accurate to a tolerance $\epsilon$ for all $t$. The resulting fast time-stepping algorithm uses this quadrature with time-independent nodes and weights, avoiding repeated re-computation. This algorithm does not impose any constraints on the final simulation time, making it suitable for modeling systems over very long durations. The storage requirement of Li’s algorithm is $O(Q)$ and its complexity is $O(NQ)$, where $N$ is the number of time steps.

Beylkin and Monzón [4] introduced an algorithmic framework based on Prony’s method and Hankel matrices for approximating the power function $t^{-\beta}$, $\beta>0$, on a finite interval by a short sum of exponentials. Their method produces, for any prescribed accuracy, a representation with far fewer terms than the standard Fourier representations would require. In a follow-up study [5], they achieved approximations for such functions with uniform relative error on the whole real line. In their method, they applied Poisson summation to discretize the integral representations of these functions. They also introduced a straightforward algorithm based on Prony’s method to remove those terms with small exponents from the initial approximations.

Building on the Beylkin-Monzón approach, McLean [14] revisited the exponential-sum approximation of $t^{-\beta}$, $\beta>0$. In order to achieve a desired accuracy, he compared the performance of the exponential-sum approximation of the Laplace integral representation using the substitutions $p=\mathrm{e}^{x}$ and $p=\exp(x-\mathrm{e}^{-x})$ both before and after applying Prony’s method. He demonstrated that the second substitution results in significantly fewer terms prior to the application of Prony’s method, while both approximations yield approximately the same number of exponential terms once Prony’s method is applied.

In this paper, we develop an efficient exponential-sum approximation for the power-law kernel $t^{\alpha-1}$ in the given interval with given fractional order $\alpha$, and $L$ exponential terms. Specifically, we pre-specify the computational cost (i.e. the number of terms $L$ in the sum of exponentials) and then minimize the approximation error. This approach differs from existing relative error tolerance-driven approaches that compute quadrature ranges and hence the number of exponential terms $L$ to meet a prescribed error and accept the resulting computational cost. Starting from the Laplace-transform integral for $t^{\alpha-1}$, we apply a composite trapezoidal rule to approximate it on $[\delta,T]$, truncating the infinite tails to render the computation finite. We then introduce Prony’s method to reduce the number of exponential terms, thereby decreasing memory usage and computational time. The result is a compact exponential-sum representation of the power-law kernel, matching or improving upon the initial maximum error before applying Prony’s method. Further, we embed this compressed kernel into Riemann-Liouville fractional integral for solving fractional differential equations, and present two solution strategies based on constant interpolation and first-order ODE formulations. The reduced ODE system is solved with the backward Euler method or the trapezoidal rule. Numerical experiments demonstrate that the proposed methods are practical and advantageous with significantly reduced computational cost for solving a wide variety of problems in fractional calculus. We also show the accuracy and effectiveness of the different solution schemes depends upon fractional order $\alpha$ and number of exponentials $L$ and discuss experimental order of convergence (EOC) for the tested schemes.

2. Kernel Approximation by Exponential Sums

2.1. The General Idea

Our aim is to construct an efficient numerical approximation by a sum of exponentials to the kernel function $t^{\alpha-1}$ of Riemann-Liouville fractional integrals of order $\alpha\in(0,1)$ for functions $f\in C[0,T]$ (where $T$ is a real positive number) with a starting point at [math] represented as

[TABLE]

The convolution kernel function $t^{\alpha-1}$ of this operator can be written as the Laplace transform of a power function in the form

[TABLE]

Choosing $p=\mathrm{e}^{\omega}$, eq. (2) reads

[TABLE]

Since the integral $(\ref{KerExp1})$ is over the unbounded domain $(-\infty,\infty)$, we truncate the domain to a finite interval $[l_{\min},l_{\max}]$, where $l_{\min}<0$ and $l_{\max}>0$ are real numbers. Then for some given $L\in\mathbb{N}$, by generating the nodes $\omega_{l}=l_{\min}+(l-1)h$ for $l=1,\ldots,L$, we divide this interval into $L-1$ subintervals of width $h=(l_{\max}-l_{\min})/(L-1)$, and finally we apply the trapezoidal rule to obtain

[TABLE]

where

[TABLE]

A finite exponential sum approximation for the kernel function can then be obtained as

[TABLE]

where $0<\delta<T$, $\theta=(w_{1},\dots,w_{L},b_{1},\dots,b_{L})\in\mathbb{R}^{2L}$ is the vector of parameters $w_{l}$ and $b_{l}$ called the weights and exponents, respectively, and

[TABLE]

With this approach, our goal is to find a set of values $\theta=(w_{1},\dots,w_{L},b_{1},\dots,b_{L})\in\mathbb{R}^{2L}$ such that the absolute error

[TABLE]

is small in modulus for $t\in[\delta,T]$, where data $\alpha,\delta,T$ and $L$ are given. As such, our course of action differs from the approach of McLean [14] who computed suitable values for $l_{\min}$ and $l_{\max}$ and hence, indirectly, $L$ based on a given relative error tolerance for the finite sum approximation $(\ref{TrapExpFin})$. Thus, while McLean intends to achieve a certain given accuracy with whatever computational cost is necessary, we pre-specify the computational cost that we are prepared to invest (which essentially depends on $L$) and want to minimize the error.

Remark 1*.*

The case of interest for us here is when $\alpha\in(0,1)$. Therefore, the function $t^{\alpha-1}$ is unbounded as $t\to 0$, a feature that the approximation $f(t,\theta)$ cannot reproduce for any $\theta$. This is the reason why, as indicated above, we aim for an approximation of $t^{\alpha-1}$ only on the interval $[\delta,T]$ with some $\delta>0$, not on the entire interval $(0,T]$. We will explain in Section 4 how we deal with the interval $(0,\delta]$ that is not covered by this approach.

Remark 2*.*

Note also that we can apply a simple rescaling technique to transfer our results from a general interval $[\delta,T]$ to the interval $[\delta/T,1]$ which is normalized in the sense that its right end point is fixed at $1$ (or vice versa): If $\theta=(w_{1},\dots,w_{L},b_{1},\dots,b_{L})$ is a choice of the weights and exponents for the interval $[\delta,T]$ and

[TABLE]

then $\tilde{\theta}=(\tilde{w}_{1},\dots,\tilde{w}_{L},\tilde{b}_{1},\dots,\tilde{b}_{L})$ can be used on $[\delta/T,1]$, and the associated errors are connected to each other via

[TABLE]

where $t\in[\delta,T]$ and $\tilde{t}=t/T\in[\delta/T,1]$.

The first technical step for us is to determine $l_{\min}$ and $l_{\max}$ for truncating the integral $(\ref{KerExp1})$ such that the integral’s contribution outside the interval $[l_{\min},l_{\max}]$ is negligible. We need the integrals

[TABLE]

to be less than a small threshold $\epsilon$.

2.2. Condition for $I_{\min}\leq\epsilon$

To find a suitable choice for $l_{\min}$, we consider the integral $I_{\min}$ in $(\ref{3})$. By substituting $u=\mathrm{e}^{\omega}$, and $v=tu$, we get

[TABLE]

where $\gamma$ is the lower incomplete Gamma function defined as

[TABLE]

see [15, §8.2]. For small $s=t\mathrm{e}^{l_{\min}}$, the function $\gamma$ has the asymptotic behavior $\gamma(1-\alpha,s)\sim s^{1-\alpha}/(1-\alpha)$ by [15, eqs. (8.2.6) and (8.7.1)]. Thus, $I_{\min}\approx\frac{\exp((1-\alpha)l_{\min})}{1-\alpha}$. Since we need $I_{\min}\leq\epsilon$, we obtain

[TABLE]

Additionally, we must make sure that $s=t\mathrm{e}^{l_{\min}}$ is sufficiently small. Hence, we must suppose that $t\mathrm{e}^{l_{\min}}<\epsilon$ for all $t\in[\delta,T]$, i.e. $T\mathrm{e}^{l_{\min}}<\epsilon$. We thus also require

[TABLE]

Combining these two inequalities, we see that we can choose a lower boundary for truncating the integral $(\ref{KerExp1})$ in the form

[TABLE]

We show in Figure 1 how different values of $\alpha$ and $T$ can affect $l_{\min}$. Figure 1(a) shows the behaviour of $f_{1}(\alpha):=\frac{\ln(\epsilon(1-\alpha))}{1-\alpha}$ for different values of $\alpha$ in $[0.01,0.99]$ and $\epsilon=10^{-8}$. Figure 1(b) shows the behaviour of $f_{2}(T):=\ln\left(\epsilon/T\right)$ for different values of $T$ in $[1,10^{4}]$ and the same $\epsilon$. We here see that different values of $T$ from the interval $[1,10^{4}]$ can have an effect on our choice for $l_{\min}$ if $\alpha$ in $[0.01,0.3]$ because both values of $f_{1}$ and $f_{2}$ change in the interval $[-18,-28]$. But for values of $\alpha$ greater than $0.3$, $T$ does not affect $l_{\min}$. Thus, the value of $\alpha$ is highly significant for approximating $l_{\min}$. It should be noted that values of $\alpha$ closer to $1$ result in more negative $l_{\min}$ and for fixed values of $\delta$ and $\epsilon$ a larger truncated interval $[l_{\min},l_{\max}]$ is then obtained.

2.3.  Condition for $I_{\max}\leq\epsilon$

Next, we consider the integral $I_{\max}$ in $(\ref{3})$. Similarly as above, by substituting $u=\mathrm{e}^{\omega}$ and $v=tu$, we get

[TABLE]

where $\Gamma$ is the upper incomplete Gamma function [15, §8.2] defined as

[TABLE]

For large $s=t\mathrm{e}^{l_{\max}}$, the function $\Gamma$ asymptotically behaves as $\Gamma(1-\alpha,s)\sim s^{-\alpha}\mathrm{e}^{-s}$, see [15, § 8.11(i)]. Thus

[TABLE]

Now, $I_{\max}\leq\epsilon$ results in $-\ln t-\alpha l_{\max}-t\mathrm{e}^{l_{\max}}\leq\ln\epsilon$. Since the dominant term is $t\mathrm{e}^{l_{\max}}$, an estimate for $l_{\max}$ dependent on $t$ is then obtained as

[TABLE]

Since for $t\geq\delta$, $l_{\max}(t)\leq l_{\max}(\delta)$, we can choose an upper bound for truncating the integral at $t=\delta$, viz.

[TABLE]

Remark 3*.*

In the case $\delta>\ln(1/\epsilon)$, we obtain $l_{\max}<0$. While this case is theoretically possible, in practical numerical settings one typically has $\delta\ll\ln(1/\epsilon)$, so such a situation is unlikely to occur. Moreover, even if it does arise, the algorithm remains well defined and may in fact yield an especially strong reduction in the number of terms after applying Prony’s method. (See also Remark 7.)

Figure 2(a) shows the behaviour of $g(\omega)$ as defined in eq. (5) for different values of $\alpha$ in $(0,1)$ at $t=1$. In all cases, the function decays very rapidly as one moves from the local maximum of $g$ to the right. Thus, we may always choose a rather small positive value for the truncation point $l_{\max}$, independent of $\alpha$. This is because the $\alpha$-dependent term in the asymptotic limit of the upper incomplete Gamma function is subdominant. Figure 2(b) shows the behaviour of $g(\omega)$ for different values of $t$ in $[\delta,T]$ and for a fixed value of $\alpha$. The larger values of $l_{\max}$ obtained by decreasing $t$ such that its maximum value is at $t=\delta$. We here choose this value, that is, $\ln\left(-\frac{\ln\epsilon}{\delta}\right)$ independent of $t$ as an approximation of $l_{\max}$. So, the value of $\delta$ plays a significant role in determining $l_{\max}$. It should be noted that by decreasing $\delta$, $l_{\max}$ become more positive and for fixed values of $\alpha$ and $\epsilon$ a larger truncated interval $[l_{\min},l_{\max}]$ is then obtained.

Remark 4*.*

As another approach, one can make an estimate of $l_{\min}$ and $l_{\max}$ by analyzing the integrand $g(\omega)=\mathrm{e}^{(1-\alpha)\omega}\mathrm{e}^{-t\mathrm{e}^{\omega}}$ to ensure that its contribution outside $[l_{\min},l_{\max}]$ is negligible, i.e., we require $\mathrm{e}^{(1-\alpha)\omega-t\mathrm{e}^{\omega}}<\epsilon$. For large positive $\omega$, $\mathrm{e}^{\omega}$ grows rapidly, and the term $-t\mathrm{e}^{\omega}$ dominates. Thus $\omega>\ln\left(-\frac{\ln\epsilon}{\delta}\right)$ and this yields an approximation of $l_{\max}$ as it is given by $(\ref{lmin_approx2})$. For large negative $\omega$, $\mathrm{e}^{\omega}\to 0$ leads to $g(\omega)\approx\mathrm{e}^{(1-\alpha)\omega}$, and then we get $\omega<\frac{\ln\epsilon}{1-\alpha}$. Thus, we obtain an approximation as $l_{\min}\approx\frac{\ln\epsilon}{1-\alpha}$. This estimate for $l_{\min}$ does not depend on $T$ but shows that it is effectively dependent on $\alpha$ similar to the first approach. It should be noted that the estimate $\frac{\ln\epsilon}{1-\alpha}$ for different values of $\alpha$ and a fixed value of $\epsilon$ exhibits similar behavior to that shown in Figure 1(a), although with very small differences in their values. However, in the numerical examples presented here, we determine the values of $l_{\min}$ and $l_{\max}$ using the first approach.

3. Prony’s Method

To improve the level of accuracy in the kernel approximation, a straightforward idea is to increase the value of $L$, i.e. to use more terms in the sum of exponentials. However, using larger values of $L$ leads to more evaluations and hence an increase in memory consumption and computation time. Thus, it is essential to equip our numerical algorithm for approximating the kernel with efficient methods that can reduce the number of exponentials in the sum while preserving accuracy.

Originally introduced in [18] and discussed in general in [19] and with an emphasis on our specific application in [4, 5, 14], the Prony’s method allows us to reduce the number of terms in the sum on the right-hand side of eq. (6), i.e. our approximation of the kernel of the integral operator, without losing accuracy. Specifically, we split up the sum as

[TABLE]

and replace the first sum on the right-hand side by an expression with the same structure but a smaller number of summands while not changing the second sum. This results in an approximation of the kernel function using fewer exponential terms. Figure 3 illustrates the behavior of weights and exponents $(w_{l},b_{l})$ for $l=1,\ldots,256$ versus the nodes $\omega_{l}$ on the truncation domain $[l_{\min},l_{\max}]$ obtained by approximating $t^{\alpha-1}$ over the interval $[10^{-2},1]$ for various values of $\alpha$ with $\epsilon=10^{-10}$. With $\omega_{l}$ in the interval $[l_{\min},0]$ the values of $b_{l}$ for $\alpha=0.1$ range from $-6.8873\times 10^{-12}$ to $-0.8806$; for $\alpha=0.5$, from $-2.5\times 10^{-21}$ to $-0.9528$; and for $\alpha=0.9$, from $-9.9\times 10^{-111}$ to $-0.6394$. As we see, the exponents $b_{l}=-\mathrm{e}^{\omega_{l}}$, where $\omega_{l}$ is within the interval $[l_{\min},0]$, mostly take values that are very near to zero in the interval $(-1,0)$. This allows us to justify the application of Prony’s method, resulting in an effective reduction. Furthermore, the weights begin with very small positive values and increase, and their maximum values decrease as $\alpha$ increases. For example,

[TABLE]

For a precise description of Prony’s method, we follow McLean [14]. Recalling first that $\omega_{l}=l_{\min}+(l-1)(l_{\max}-l_{\min})/(L-1)$ for $l=1,2,\ldots,L$, we set

[TABLE]

so that $M\leq L$, $\omega_{l}\leq 0$ for $l=1,2,\ldots,M$ and $\omega_{M+1}>0$ if $M<L$. Furthermore, we choose some $L_{\text{p}}\in\{1,2,\ldots,M\}$ and a positive integer $K$ such that $2K-1\leq L_{\text{p}}$. The goal of Prony’s method is to replace the first $L_{\text{p}}$ summands in the right-hand side of eq. (6) by

[TABLE]

i.e. by a sum with at most a bit more than half as many (and, if possible, even fewer) terms as the original one, where the parameters $\rho_{k}$ and $\eta_{k}$ are chosen such that we can maintain the desired accuracy, i.e. such that

[TABLE]

is small enough in modulus for all $t\in[\delta,T]$, thus permitting to reduce the total number of terms from $L$ to

[TABLE]

Based on Prony’s method, in order to find the $2K$ parameters $\rho_{k}$ and $\eta_{k}$, we take the following steps:

Check that $2K-1\leq L_{\text{p}}$ and compute

[TABLE] 2. 2.

Construct the $K\times K$ Hankel matrix

[TABLE]

and the vector $\mathbf{b}=-(g_{K},g_{K+1},\ldots,g_{2K-1})^{\mathrm{T}}$, and solve the linear system $\mathbf{H}\mathbf{q}=\mathbf{b}$ with $\mathbf{q}=(q_{0},q_{1},\ldots,q_{K-1})^{\mathrm{T}}$. 3. 3.

Form the monic polynomial $q(z)=q_{0}+q_{1}z+\ldots+q_{K-1}z^{K-1}+z^{K}$, and find its roots $\eta_{k}$ ($k=1,\ldots,K$) to obtain the Prony’s exponents for eq. (19). 4. 4.

Construct the $(2K)\times K$ Vandermonde matrix

[TABLE]

and solve the linear system $\mathbf{V}\boldsymbol{\rho}=\mathbf{g}$ where $\mathbf{g}=(g_{0},g_{1},\ldots,g_{2K-1})^{\mathrm{T}}$ and $\boldsymbol{\rho}=(\rho_{1},\rho_{2},\ldots,\rho_{k})^{\mathrm{T}}$ to determine the Prony’s weights $\rho_{k}$ in eq. (19).

Remark 5*.*

The roots of the polynomial $q(z)$ that arises in item 3 above might not all be real and positive. In addition, the $K\times K$ Hankel matrix mentioned in item 2 might be nearly singular in some cases. In practice, the detailed discussion in [4, 5, 14] shows that these issues only cause minor difficulties for our algorithm that can easily be circumvented.

Remark 6*.*

Although only the first $K$ rows of the Vandermonde system in item 4 would suffice to uniquely determine the weights $\rho_{k}$, in practice this choice may lead to severe numerical ill-conditioning, especially when the nodes $\eta_{k}$ are very small. By formulating the problem as an overdetermined $(2K)\times K$ system, we gain additional stability and robustness against numerical errors, while still recovering the same solution in well-conditioned cases.

To determine suitable values for $K$ and $L_{\text{p}}$, we employ the method outlined in [14, Example 1] but present it as an algorithm. Assume that the original approximation $t^{\alpha-1}\approx\sum_{l=1}^{L}w_{l}\mathrm{e}^{b_{l}t}$ yields a maximum error of $\epsilon^{\prime}$ over the interval $[\delta,T]$. The steps of the method are as follows:

Initialize with $K=1$ and $L_{\text{p}}=M$ to maximize the initial reduction in the number of points. 2. 2.

Compute the parameters of Prony’s method:

• •

Calculate the exponents $\eta_{k}$ and weights $\rho_{k}$ for $k=1,\dots,K$ using Prony’s method as described above.

• •

Evaluate the Prony approximation error $\varepsilon_{\text{p}}(t)$ from $(\ref{pronyerror})$ for $t\in[\delta,T]$. 3. 3.

Check the error condition:

• •

If $\max_{t\in[\delta,T]}\varepsilon_{\text{p}}(t)\leq\epsilon^{\prime}$, accept the current values of $K$ and $L_{\text{p}}$ as suitable choices for reducing the total number of terms. Terminate the algorithm.

• •

Otherwise, proceed to step $4$. 4. 4.

Update parameters:

• •

If $2K-1<L_{\text{p}}$, decrement $L_{\text{p}}$ by 1 and go back to step 2.

• •

If $2K-1\geq L_{\text{p}}$, increment $K$ by 1, reset $L_{\text{p}}=M$, and go back to step 2, provided $2K-1<L_{\text{p}}$. Repeat until the error condition is satisfied.

With these considerations, we can replace the initial sum-of-exponentials approximation $f(t,\theta_{0})$, say, for $t^{\alpha-1}$ from eq. (6) with a certain parameter vector $\theta_{0}$ by the new, computationally cheaper, approximation

[TABLE]

where $\theta_{\text{p}}=(\rho_{1},\rho_{2},\ldots,\rho_{K},w_{L_{\text{p}}+1},w_{L_{\text{p}}+2},\ldots,w_{L},\eta_{1},\eta_{2},\ldots,\eta_{K},b_{L_{\text{p}}+1},b_{L_{\text{p}}+2},\ldots,b_{L})$. For this final approximation after the Prony reduction step, we can immediately derive an upper bound for the absolute value of the error $e_{\text{p}}(t)=t^{\alpha-1}-f_{\text{p}}(t,\theta_{\text{p}})$ according to

[TABLE]

This gives us $\max_{t\in[\delta,T]}|e_{\text{p}}(t)|\leq\epsilon^{\prime}+\epsilon^{\prime}=2\epsilon^{\prime}$. Therefore, after reduction, this algorithm allows us to preserve the maximal error close to the initial maximal error obtained before applying Prony’s method.

Table 1 presents a comparison of the maximum error values $\max_{t\in[10^{-2},1]}|e(t,\theta_{0})|$ before applying Prony’s method with the corresponding values $\max_{t\in[10^{-2},1]}|e_{\text{p}}(t)|$ after the Prony procedure for various values of $L$ from $32$ to $1024$ and $\alpha\in\{0.1,0.5,0.9\}$. The results in Table 1 show the following:

In most cases, the maximum error after applying Prony’s method, viz. $\max_{t\in[10^{-2},1]}|e_{\text{p}}(t)|$, either does not change or actually improves upon the initial error $\max_{t\in[10^{-2},1]}|e(t,\theta_{0})|$. For example,

• •

in the case $\alpha=0.1$ and $L=32$, the error slightly decreases from $5.102467\times 10^{-2}$ to $5.102386\times 10^{-2}$,

• •

for $\alpha=0.1$ and $L=128$, the error significantly improves from $1.320726\times 10^{-8}$ to $1.980379\times 10^{-10}$, and

• •

for $\alpha=0.5$ and $L=256$, the error decreases from $3.518998\times 10^{-10}$ to $5.593037\times 10^{-11}$.

This confirms that the algorithm maintains the accuracy within the theoretical bound $2\epsilon^{\prime}$ or even enhances it. 2. 2.

The value of $K$ is always small (ranging from 1 to 5). This indicates an efficient compression of the original exponential sum. For instance, for $\alpha=0.9,L=128$, only $K=2$ terms are needed, reducing $L_{\text{f}}$ to 6, demonstrating that a minimal number of terms can achieve comparable accuracy. 3. 3.

The algorithm initializes $L_{\text{p}}=M$, and the iterative process adjusts $L_{\text{p}}$ and $K$ to satisfy the error criterion $\max_{t\in[10^{-2},1]}|\varepsilon_{p}(t)|\leq\epsilon^{\prime}$. The table shows that in all cases, the final values of $L_{\text{p}}$ are equal to $M$, indicating that the algorithm can use all available negative nodes initially, optimizing the reduction process. 4. 4.

For all $\alpha$, the algorithm significantly reduces the number of terms ($L_{\text{f}}\ll L$). As $\alpha$ increases, $L_{\text{f}}$ becomes particularly small, i.e. especially few terms are required for comparable accuracy. Thus, Prony’s method effectively reduces the number of terms in the exponential sum approximation.

Therefore, Table 1 confirms that the proposed algorithm effectively reduces the number of terms in the exponential sum approximation while keeping the maximum error close to or better than the initial error. The small values of $K$ and the significant reduction of the total number of terms, and hence the computational cost, from $L$ to $L_{\text{f}}$ highlight the algorithm’s efficiency, making it highly suitable for applications where computational efficiency and accuracy are critical.

In all tables, the notation $5.102467(-2)$ with a mantissa and an exponent in parentheses represents the number $5.102467\times 10^{-2}$.

Remark 7*.*

If the parameters introduced in Section 2 satisfy $\delta\geq\ln(1/\epsilon)$ then, by eq. (17), we have $\omega_{L}=l_{\max}\leq 0$ and hence eq. (18) implies $M=L$. It is not uncommon for the algorithm described above (see the examples shown in Table 1) to yield $L_{\text{p}}=M=L$, so that in this case the second sum on the right-hand side of eq. (21) is empty, thus we obtain an especially large reduction in the number of terms.

For a more precise demonstration of the performance of the Prony method, we visualize in Figure 4 a comparison of the difference between the exact and approximate values of the power law kernel before the application of Prony’s method, i.e. $e(t,\theta_{0})$, vs. $t$ (plotted as continuous lines) with the analog error after using Prony’s method, viz. $e_{\text{p}}(t)$, vs. $t$ (dashed) for various values of $L$ and $\alpha\in\{0.1,0.5,0.9\}$, where $t$ ranges over the interval $[10^{-2},1]$. As we see, for smaller values of $L$, with Prony’s method the difference grows as $t$ increases, but the maximum absolute error remains constant. As $L$ increases, the difference exhibits the same behaviour as observed before Prony’s method.

Remark 8*.*

To ensure stability in Prony’s reduction algorithm, the Hankel matrix $\mathbf{H}$ must be sufficiently far away from a singular matrix. Based on our numerical observations, for $T=1$, this condition holds, making the algorithm stable. However, as $T$ increases, the Hankel matrix may become singular in some cases, causing algorithm instability. To address this situation in Prony’s method, the rescaling approach discussed in Remark 2 can be applied. Therefore, we first compute the exponents and weights for the interval $[\delta/T,1]$, then apply Prony’s algorithm, and finally rescale as in eq. (9) to create an approximation for the interval $[\delta,T]$ with a reduced number of terms. For more details on the maximum error values in this case, see Appendix.

In Algorithm 1, we now present our algorithm of Prony’s method formally in a pseudo-code style.

The algorithm inputs include data $\alpha$ (fractional order), $\delta$ (initial time), $T$ (final time), $N$ (number of nodes for the discretization of the time interval $[\delta,T]$), $L$ (number of exponential terms), and $\epsilon$ (a very small threshold for neglecting the integrals (11)). This algorithm is split into two parts:

• •

Part 1 (Initialization)—lines 1–20: initializes the input data, then computes the lower and upper bounds of the integral $(\ref{KerExp1})$ based on equations $(\ref{8})$ and $(\ref{lmin_approx2})$ and initial values of the weights and exponents $\theta_{0}=[w_{0},b_{0}]\in\mathbb{R}^{2L}$, where $w_{0}=[w_{1}^{0},\dots,w_{L}^{0}]$ and $b_{0}=[b_{1}^{0},\dots,b_{L}^{0}]\in\mathbb{R}^{L}$ by $(\ref{weightandnode})$, and finally evaluates the initial approximation error of the kernel function, $e(t,\theta_{0})$ using $(\ref{Abserror})$ and its maximum error $\epsilon^{\prime}$.

• •

Part 2 (Application of Prony’s method)—lines 22–31: applies the Prony’s method to the initial approximation $\theta_{0}$ to reduce the total number of exponential terms while preserving the accuracy $\epsilon^{\prime}$.

Remark 9*.*

One can choose a uniform discretization of the interval $[\delta,T]$ on which the approximation errors are computed. However, to reduce the computation time, especially for large $T$, we can use a geometric grid in the interval $[\delta,T]$ given by $t_{j}=\delta(T/\delta)^{(j-1)/(N-1)}$ for $j=1,\ldots,N$. In that case, the common ratio between consecutive points is $r=\left(T/\delta\right)^{1/(N-1)}$.

4. Application to Fractional Ordinary Differential Equations

Next, our aim is to present two algorithms for the numerical solution of the fractional ordinary differential equation (FODE)

[TABLE]

subject to the initial condition $y(0)=y_{0}$, where $0<\alpha<1$, $T$ is a real positive number, and $D^{\alpha}$ is the Caputo type fractional derivative [10, Chapter 3] of order $\alpha$, which is defined by

[TABLE]

where $I^{1-\alpha}$ is the Riemann-Liouville (RL) fractional integral of order $1-\alpha$ defined by $(\ref{RLIntegral})$ which leads to

[TABLE]

The Caputo type fractional derivative $D^{\alpha}$ of order $\alpha$ is a left inverse of the RL integral $I^{\alpha}$ of order $\alpha$, that is, $D^{\alpha}I^{\alpha}y=y$, but it is not a right inverse since

[TABLE]

By applying the RL fractional integral of order $\alpha$ to both sides of Equation (23), and considering Equation (26), we get the Volterra integral equation

[TABLE]

see [10, § 6.1]. To numerically solve (27), we first introduce an arbitrary discretization (e.g., a uniform grid), which is defined as $0=t_{0}<t_{1}<\ldots<t_{N}=T$, where the step sizes are $h_{j}=t_{j}-t_{j-1}$ for $j=1,\ldots,N$, and $h_{j}$ may vary. Here, we assume that the approximations $y_{j}\approx y(t_{j})$ for $j=1,2,\ldots,n-1$ and some $n\geq 1$ have been already computed. Then, we compute the approximation $y_{n}\approx y(t_{n})$ using the equation

[TABLE]

where $\tilde{I}^{\alpha}f(t_{n},y_{n})$ is a suitable approximation for the RL integral

[TABLE]

computed with the help of the algorithm described in Sections 2 and 3. Following the approach outlined in [13, 9], we split up the integral in equation (29) into a local part $L^{\alpha}_{n}$ and a history part $H^{\alpha}_{n}$ according to

[TABLE]

where the history part is given by

[TABLE]

and the local part is defined by

[TABLE]

Here, to approximate (31) and (32), we present two algorithms.

4.1. Method of Constant Interpolation

In the first approach, to approximate the integrand $f$ in (31) and (32), we use piecewise constant interpolation $f(\tau,y(\tau))\approx f(t_{j},y_{j})$ for all $\tau\in(t_{j-1},t_{j})$ and all $j$.

In the local part $L^{\alpha}_{n}$ the integration ranges from $t_{n-1}$ to $t_{n}$, i.e., it does not involve contributions from the process history. Due to the singularity of the kernel function $(t_{n}-\tau)^{\alpha-1}$ at $\tau=t_{n}$, there is a region near $\tau\approx t_{n}$ where the exponential sum approximation loses accuracy. Therefore, to compute the local contribution, we need to use standard numerical methods that avoid relying on the exponential sum approximation, achieving accurate results with $O(1)$ time complexity and $O(1)$ memory usage. Thus, we compute $L^{\alpha}f(t_{n},y_{n})$ as

[TABLE]

where

[TABLE]

Here, we have used the functional equation $\alpha\Gamma(\alpha)=\Gamma(\alpha+1)$ [1].

In the history part $H^{\alpha}_{n}$, the integration ranges from [math] to $t_{n-1}$, meaning the history part is confined to times up to at most $t_{n-1}$ and thus it (a) avoids the kernel function’s singularity at $\tau=t_{n}$, and (b) becomes zero for $n=1$. As a result, the exponential sum technique can be used to efficiently compute an approximation for $H^{\alpha}_{n}[f(\cdot,y(\cdot))]$. By considering the time discretization, we obtain

[TABLE]

Using the piecewise constant interpolant $f(\tau,y(\tau))\approx f(t_{j},y_{j})$ over each interval $[t_{j-1},t_{j}]$ for $j=1,\dots,n-1$, and considering the approximation of the kernel function $(t_{n}-\tau)^{\alpha-1}$ for $0\leq\tau\leq t_{n-1}$ as given in (6) and applying it to the sum, we get

[TABLE]

where $c_{\alpha}=\frac{1}{\Gamma(\alpha)\Gamma(1-\alpha)}=\frac{\sin\pi\alpha}{\pi}$ in view of the reflection formula (one of the fundamental properties of the Gamma function [1]) and

[TABLE]

with

[TABLE]

It should be noted that the value of $\delta$ used in the construction of the sum-of-exponentials approximation for the kernel function must be less than or equal to the smallest time step $\min_{j}h_{j}=\min_{j}(t_{j}-t_{j-1})$. This ensures that $\delta\leq h_{n}=t_{n}-t_{n-1}\leq t_{n}-\tau\leq T$ holds for all $0\leq\tau\leq t_{n-1}$, which implies $\delta\leq t_{n}-\tau\leq T$.

The crucial property of this approach is that we can find a relation between $K_{l,n,j}$ and $K_{l,n-1,j}$, namely

[TABLE]

Thus, in this evaluation scheme, we derive the recursive relation

[TABLE]

with the initialization $\Phi_{l}^{1}=0$. Using eq. (38) instead of eq. (35), we can reduce the computational complexity of the evaluation of $\Phi_{l}^{n}$ (required in the $n$th time step) from $O(n)$ to $O(1)$, and therefore the overall cost for taking steps $n=1,2,\ldots,N$ goes down from $O(N^{2})$ to $O(N)$. The approximate solution to the fractional differential equation at $t_{n}$ is thus computed as

[TABLE]

Here, we estimate the value of $y(t)$ at time $t_{n}$ using only the information from the immediately preceding time point $t_{n-1}$, without relying on any data from the earlier history. Thus, this approach can also efficiently reduce the active storage requirements because we do not need to keep any history data in memory for the later time steps..

Next we need to solve eq. (39) for the unknown variable $y_{n}$ that appears on both sides. But the nonlinearity of the function $f$ makes it generally impossible to find $y_{n}$ directly. Thus, we have to resort to an iterative method for solving this equation. Here we use the iterative Newton-Raphson method, which starts with an initial estimate $y_{n}^{(0)}$ for $y_{n}$ and then computes the successive approximations of $y_{n}$ using the following equation

[TABLE]

where

[TABLE]

Here, $y_{n}^{(s)}$ is the $s$-th approximation of $y_{n}$, with the initial guess $y_{n}^{(0)}=y_{n-1}$. This is a common approach—it is often effective to begin Newton iterations with the most recently calculated approximation. We then terminate our iterations when two successive values differ by less than a specified tolerance.

Remark 10*.*

If the Newton iteration does not converge (which may happen, e.g., if the initial guess $y_{n}^{(0)}$ is not good enough), or if the computation of the function $G^{\prime}$ is infeasible or too expensive, then one can also use a fixed point iteration directly based on eq. (39), i.e. one (as above) starts with $y_{n}^{(0)}=y_{n-1}$ and then proceeds according to

[TABLE]

again terminating when two successive values are sufficiently close to each other and then accepting the last one as $y_{n}$. Typically, this procedure will require more iterations to converge than the Newton-Raphson method, but each step is computationally cheaper.

It should be noted that we can easily compute $K_{l,n,n-1}$ formally using (36) which yields

[TABLE]

As noted in Section 3, before applying Prony’s method, we are likely to have exponents $b_{l}$ near zero, which causes a cancellation of significant digits or even a division by zero when using (41) naively and thus may be fatal when implementing the algorithm in finite-precision arithmetic. However, Prony’s method can generate negative exponents, which are not very close to zero, particularly for smaller intervals, enabling the use of equation (41) without causing any problems, and for very large intervals, where very small exponents may exist, exact values of $K_{l,n,n-1}$ by (41) can then be replaced by numerical computation of $K_{l,n,n-1}$ using integration methods (see Table 4 in Appendix).

4.2. Method of First Order ODE

Here, for the local part, we use a linear interpolant as the straight line between the points $(t_{n},f_{n})$ and $(t_{n-1},f_{n-1})$. So, we get the approximation $f(\tau,y(\tau))\approx(f_{n}-f_{n-1})(\tau-t_{n-1})/h_{n}+f_{n-1}$ for all $\tau\in(t_{n},t_{n-1})$. By substituting this into (32) and after some calculation, we get

[TABLE]

For the history part, after considering the approximation of the kernel function, we have

[TABLE]

where

[TABLE]

It is easy to show that for any $b_{l}$ with $l=1,\ldots,L$, the function $\mu(\cdot,h,b_{l})$ is the unique solution to the first order differential equation

[TABLE]

subject to the initial condition $\mu(h,h,b_{l})=0$.

Now, we require a method for the approximate solution to the initial value problem (45). Since all the $b_{l}$ are negative and some of them have large absolute values, i.e. $b_{l}\ll 0$, the differential equation (45) becomes very stiff, posing challenges for numerical solutions. Therefore, an A-stable scheme is necessary, and preferably one should even choose an L-stable method. Implicit methods are preferred to achieve reliable results [17]. Here, we will test both the backward Euler method (which is L-stable) and the trapezoidal rule (that is not L-stable) for this purpose. Notably, equation (45) indicates that the initial value problem must be solved numerically $L$ times, corresponding to each parameter value $b_{l}$ ($l=1,2,\dots,L$).

Thus, for $n\geq 1$, we get the backward Euler scheme, which is of first order

[TABLE]

and the trapezoidal rule, which is of second order

[TABLE]

with $\mu(t_{1},h_{1},b_{l})=\mu(h_{1},h_{1},b_{l})=0$. Thus, the approximate solution to the fractional differential equation (23) at $t_{n}$ is

[TABLE]

Similar to the method of constant interpolation, solving (48) for $y_{n}$ is here done by Newton’s iteration method or a direct fixed point iteration.

5. Numerical Results

Here we present some numerical examples of fractional differential equations and solve them using the methods introduced in Section 4 and our new optimization algorithm for approximating the kernel function by exponential sums. Our algorithm has been implemented in MATLAB R2023b on a notebook with an AMD Ryzen 7 PRO 7840U CPU clocked at 3.30 GHz running Windows 11. In all examples, the values of $\epsilon$ and the tolerance for the Newton iterations are set to $10^{-10}$ and for the kernel approximation, $\delta$ is set to $10^{-5}$. For Examples 1 and 2, we assume that the discretization of the interval $[0,T]$ is uniform. i.e., $t_{n}=nh,$ where $h=T/N$ is the step size and $N\in\mathbb{N}$.

Example 1*.*

Our first example for testing our numerical method is the nonlinear problem

[TABLE]

with the initial condition $y(0)=0$ originally proposed in [11]. The exact solution of this equation is $y(t)=t^{8}-3t^{4+\frac{\alpha}{2}}+\frac{9}{4}t^{\alpha}$.

Example 2*.*

For our second example, we consider the following fractional differential equation

[TABLE]

whose exact solution is $y(t)=E_{\alpha}(\lambda t^{\alpha})$, with $E_{\alpha}$ being the Mittag-Leffler function of order $\alpha$, which can be numerically evaluated by an algorithm presented in [12]. Here, we numerically solved the problem in the interval $[0,10]$ with various values of $\alpha$ and $\lambda=-1$.

For these examples, we compare the performance of the three numerical methods (Constant Interpolation (CI), Backward Euler (BE), and Trapezoidal Rule (TR)) for solving the fractional differential equations (49) and (50), respectively, with fractional orders $\alpha\in\{0.1,0.5,0.9\}$. Tables 2–7 show the absolute errors at the final time $T$ for varying time step sizes $h=2^{-2},2^{-3},\ldots,2^{-10}$, different numbers of exponential terms $L$, and corresponding $L_{\text{f}}$ values obtained after using Prony’s method. In these tables, the column marked ‘EOC’ states the experimentally determined order of convergence for each value of $L$ which is calculated using the error at the final time by $\log_{2}(E(h)/E(h/2))$, where $E(h)=|y(T)-y_{N}|$ is the absolute error corresponding to the step size $h$, and the symbol ‘***’ indicates that EOC is not computed because the errors exhibit unexpected behavior, i.e., increasing instead of decreasing. The results of these tables indicate the following facts:

For all methods and all values of $\alpha$, absolute errors decrease as the time step size $h$ decreases (i.e., from $2^{-2}$ to $2^{-10}$). This is expected, as smaller time steps typically improve numerical accuracy. 2. 2.

The TR method consistently achieves the lowest errors and highest EOC compared to CI and BE methods for sufficiently large $L$. In Example 1, for all values of $\alpha$ and sufficiently large $L$ ($128$ for $\alpha=0.1$ and $0.5$, and $512$ for $\alpha=0.9$), EOC well matches the full second order convergence as we expect for the trapezoidal rule. In Example 2, for the TR method, the EOC is always close to $1+\alpha$. But here, increasing $L$ has not improved the EOC, see columns with $L=256$ and $1024$. Figures 5 and 6 further illustrate the high accuracy of TR compared to CI and BE methods over the entire interval for all values of $\alpha$. 3. 3.

The BE and CI methods generally exhibit values of errors and EOC that are relatively close to each other, particularly for larger values of $\alpha$, as shown in Figures 5 and 6. However, both CI and BE methods are limited to an EOC of approximately $1$ for sufficiently large $L$, indicating first-order convergence. This matches the well known fact the backward Euler is of first order. 4. 4.

Increasing $L$ often reduces errors and improves EOC for all methods and $\alpha$, indicating that more exponential terms enhance accuracy and convergence. Thus, the choice of $L$ is critical to achieve optimal convergence. However, beyond a certain $L$, further increases do not significantly reduce the absolute errors of the solution anymore. This demonstrates the efficiency of our schemes which require only a modest number of exponential terms to achieve optimal convergence. This is evident in Figures 7 and 8, where for $\alpha=0.5$ and all methods, no accuracy improvement occurs when increasing $L$ from $128$ to $256$.

Example 3*.*

Here we consider the differential equation of motion of a massless one-dimensional fractional Kelvin-Voigt model subject to the external force $f(t)=1$ for $t\neq 0$,

[TABLE]

with $c=100,k=10$ and $\alpha=0.3$, see [20]. We use the same discretization of the variable $t$ used in [20] by the mesh points

[TABLE]

with $h_{1}=10^{-4}$ and $h_{j}=1.005h_{j-1}$ for $j=2,\dots,5000$. The exact solution of this equation is

[TABLE]

where $E_{\alpha}$ is the Mittag-Leffler function of order $\alpha$. The numerical methods CI, BE, and TR are applied to solve problem (51). A comparison of absolute errors over time, as shown in Figure 9(a), reveals that the BE and TR methods exhibit similar behavior but do not significantly improve compared to the CI method. These findings align with the results from previous examples, which suggest that the TR method can offer considerable improvements over the CI and BE methods mainly for larger $\alpha$ values. Nevertheless, the trapezoidal rule (TR) can be considered as a more efficient option compared to the other two methods for solving fractional differential equations. To verify the impact of the parameter $L$ on solution accuracy, Figure 9(b) illustrates that for the trapezoidal rule, increasing $L$ enhances accuracy, but beyond $L\geq 64$, further increases yield no significant improvement. Figure 9(c) compares the numerical solutions using the trapezoidal rule for $L=16$ and $L=64$ with the exact solution from (53), demonstrating that a considerable error reduction requires a sufficiently large $L$.

6. Conclusions

In this work, we proposed a sum of exponential approximation for power functions $t^{\alpha-1},0<\alpha<1$ with pre-specified number of exponential terms. We then combined this with Prony’s algorithm to efficiently decrease the number of exponentials while preserving or even enhancing the maximum approximation error. The findings highlight the importance of choosing the appropriate method for solving fractional differential equations and number of exponential terms $L$ based on the fractional order $\alpha$ to achieve optimal accuracy and convergence. Here, the Trapezoidal Rule (TR) is the most effective method for solving fractional differential equations, offering the lowest errors and highest EOC. The CI and BE methods exhibit slower convergence and lower accuracy comparatively, rendering them suboptimal for high-precision applications. In all methods, increasing the number of exponential terms $L$ enhances the solution accuracy.

Appendix

Table 8 presents the same parameters as Table 1, but evaluated over the interval $[10^{-2},10^{3}]$. It includes the initial maximum error, $\max_{t\in[10^{-2},10^{3}]}|e(t,\theta_{0})|$, and the maximum error after applying Prony’s method, $\max_{t\in[10^{-2},10^{3}]}|e_{\text{p}}(t)|$, both of which are computed using rescaling of the weights and exponents according to Remark 2 for the interval $[10^{-5},1]$. Table 8 confirms that Prony’s method, with the rescaling approach, effectively approximates the kernel function over the large interval $[10^{-2},10^{3}]$ without loss in accuracy and produces a significant term reduction ($L_{\text{f}}\ll L$) across all cases. As in Table 1, the parameter $K$ ranges from $1$ to $5$, and in all cases, the obtained $L_{\text{p}}$ is equal to the initialized value $M$. Compared to Table 1, the value of $L_{\text{f}}$ increases for fixed $\alpha$ and $L$ because reducing $\delta$ from $10^{-2}$ to $10^{-5}$ widens the interval, resulting in a greater number of positive nodes.

Therefore, Table 8 demonstrates the robustness of Prony’s method for larger intervals, making it suitable for practical applications where computational efficiency and accuracy over wide ranges are critical.

In Table 9, we focus on fixed $L$ values ($L=256$ for $\alpha=0.1,0.5$, and $L=1024$ for $\alpha=0.9$) and varying the parameter $\epsilon$ from $10^{-10}$ to $10^{-13}$. The table shows that the final number of terms ($L_{\text{f}}$) decreases as $\epsilon$ becomes stricter (smaller). For $\alpha=0.1$, $L_{\text{f}}$ decreases from $96$ ($\epsilon=10^{-10}$) to $83$ ($\epsilon=10^{-13}$). For $\alpha=0.5$, $L_{\text{f}}$ decreases from $65$ to $55$, and for $\alpha=0.9$, $L_{\text{f}}$ decreases from $61$ to $51$. This reduction in $L_{\text{f}}$ demonstrates that Prony’s method effectively compresses the number of terms while maintaining accuracy, especially for stricter $\epsilon$ values. Furthermore, decreasing $\epsilon$ improves both initial and Prony errors, with Prony errors remaining close to or better than initial errors in most cases.

In Table 10, we have shown the values of the parameters $\rho_{k}$ (weights) and $\eta_{k}$ (exponents) after the application of Prony’s method in the cases of $L=256$ for $\alpha=0.1,0.5$ and $L=1024$ for $\alpha=0.9$ over three intervals, $[10^{-2},1]$, $[10^{-5},1]$, and $[10^{-2},10^{3}]$. In comparison with what we explained about Figure 3, before using Prony’s method, there are very small (very close to zero) weights $w_{l}$ and exponents $b_{l}$ obtained by the nodes $\omega_{i}$ on the interval $[l_{\min},0]$. However, with Prony’s method we get weights $\rho_{k}$ and exponents $\eta_{k}$ which are orders of magnitude larger. The weights $\rho_{k}$ are generally larger for smaller intervals $[10^{-2},1]$, $[10^{-5},1]$ and decrease significantly for the larger interval $[10^{-2},10^{3}]$ due to the rescaling described in Remark 2 ($w_{l}=\tilde{w}_{l}/T^{1-\alpha}$). As $\alpha$ increases, the weights are significantly larger (especially for smaller intervals). The exponents $\eta_{k}$ for the large interval $[10^{-2},10^{3}]$ are orders of magnitude smaller (closer to zero) than those for $[10^{-2},1]$ or $[10^{-5},1]$ due to the rescaling factor $b_{l}=\tilde{b}_{l}/T$.