Skip to content

Commit e08f2b8

Browse files
author
prmiles
committed
updated format and testing for fdc
1 parent 8fc16e7 commit e08f2b8

File tree

2 files changed

+74
-9
lines changed

2 files changed

+74
-9
lines changed

pyfod/fdc.py

+28-9
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,9 @@
22

33
import sys
44
import numpy as np
5-
from scipy.special import gamma as sc_gamma
65
import sympy as sp
7-
from GaussLegendreRiemannSum import GaussLegendreRiemannSum
6+
from scipy.special import gamma as sc_gamma
7+
from pyfod.GaussLegendreRiemannSum import GaussLegendreRiemannSum
88
from pyfod.GaussLegendre import GaussLegendre
99
from pyfod.GaussLaguerre import GaussLaguerre
1010
from pyfod.RiemannSum import RiemannSum
@@ -18,13 +18,13 @@ def fdc(f, start, finish, dt=1e-4, alpha=0.0, quadrature='GLegRS', **kwargs):
1818
Q2 = quad(start=start, finish=finish-dt, alpha=alpha, **kwargs)
1919
I2 = Q2.integrate(f=f)
2020

21-
if quadrature == 'GLag':
21+
if quadrature.lower() == 'glag':
2222
extend_precision = True
2323
else:
2424
extend_precision = False
2525

2626
if extend_precision is True:
27-
fd = (I1-I2)/(dt*sp.gamma(1-alpha))
27+
fd = float((I1-I2)/(dt*sp.gamma(1-alpha)))
2828
else:
2929
fd = (I1-I2)/(dt*sc_gamma(1 - alpha))
3030
# assemble output
@@ -33,13 +33,13 @@ def fdc(f, start, finish, dt=1e-4, alpha=0.0, quadrature='GLegRS', **kwargs):
3333

3434
def select_quadrature_method(quadrature):
3535
methods = dict(
36-
GLegRS=GaussLegendreRiemannSum,
37-
GLag=GaussLaguerre,
38-
GLeg=GaussLegendre,
39-
RS=RiemannSum
36+
glegrs=GaussLegendreRiemannSum,
37+
glag=GaussLaguerre,
38+
gleg=GaussLegendre,
39+
rs=RiemannSum
4040
)
4141
try:
42-
quad = methods[quadrature]
42+
quad = methods[quadrature.lower()]
4343
return quad
4444
except KeyError:
4545
print('Invalid quadrature method specified: {}'.format(quadrature))
@@ -57,6 +57,8 @@ def fcos(t):
5757
def fexp(t):
5858
return np.exp(2*t)
5959

60+
def fspexp(t):
61+
return sp.exp(2*t)
6062

6163
start = 0.0
6264
finish = 1.0
@@ -74,3 +76,20 @@ def fexp(t):
7476
alpha = 0.9
7577
out = fdc(f=fexp, alpha=alpha, start=start, finish=finish, dt=dt, NRS=NRS)
7678
print('D^{}[exp(2t)] = {} ({})'.format(alpha, out['fd'], 13.8153))
79+
80+
# Test Extended Precision - Gauss Laguerre Quadrature
81+
# Test alpha = 0.0
82+
alpha = 0.0
83+
out = fdc(f=fspexp, alpha=alpha, start=start, finish=finish, dt=dt,
84+
quadrature='glag')
85+
print('D^{}[exp(2t)] = {} ({})'.format(alpha, out['fd'], 7.38906))
86+
# Test alpha = 0.1
87+
alpha = 0.1
88+
out = fdc(f=fspexp, alpha=alpha, start=start, finish=finish, dt=dt,
89+
quadrature='glag')
90+
print('D^{}[exp(2t)] = {} ({})'.format(alpha, out['fd'], 7.95224))
91+
# Test alpha = 0.9
92+
alpha = 0.9
93+
out = fdc(f=fspexp, alpha=alpha, start=start, finish=finish, dt=dt,
94+
quadrature='glag', N=32, n_digits=60)
95+
print('D^{}[exp(2t)] = {} ({})'.format(alpha, out['fd'], 13.8153))

test/test_fdc.py

+46
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
# -*- coding: utf-8 -*-
2+
import unittest
3+
from pyfod import fdc
4+
import numpy as np
5+
import sympy as sp
6+
7+
8+
def fexp(t):
9+
return np.exp(2*t)
10+
11+
12+
def fsp(t):
13+
return sp.exp(2*t)
14+
15+
16+
# --------------------------
17+
class SelectQuadratureMethodTesting(unittest.TestCase):
18+
19+
def test_selection(self):
20+
with self.assertRaises(SystemExit):
21+
fdc.select_quadrature_method(quadrature='hello')
22+
23+
24+
# --------------------------
25+
class FDC(unittest.TestCase):
26+
27+
def check_contents(self, out):
28+
keys = ['fd', 'I1', 'I2', 'Q1', 'Q2']
29+
for key in keys:
30+
self.assertTrue(key in out,
31+
msg=str('{} not in output'.format(key)))
32+
self.assertTrue(isinstance(out['fd'], float),
33+
msg='Expect float')
34+
35+
def test_fdc(self):
36+
start = 0.0
37+
finish = 1.0
38+
NRS = 10
39+
# Test alpha = 0.0
40+
alpha = 0.0
41+
out = fdc.fdc(f=fexp, alpha=alpha, start=start,
42+
finish=finish, NRS=NRS)
43+
self.check_contents(out)
44+
out = fdc.fdc(f=fsp, alpha=alpha, start=start,
45+
finish=finish, quadrature='glag')
46+
self.check_contents(out)

0 commit comments

Comments
 (0)