Skip to content
Snippets Groups Projects
Commit 3e4f7a25 authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

add get_fornberg_coeffs

parent dc32d0dc
No related branches found
No related tags found
No related merge requests found
...@@ -71,36 +71,39 @@ def padd_with_zeros( ...@@ -71,36 +71,39 @@ def padd_with_zeros(
try: try:
import sympy as sp import sympy as sp
rational0 = sp.Rational(0)
rational1 = sp.Rational(1)
except ImportError:
rational0 = 0.0
rational1 = 1.0
def get_fornberg_coeffs( def get_fornberg_coeffs(
x = None, x = None,
a = None): a = None):
N = len(a) - 1 N = len(a) - 1
d = [] d = []
for m in range(N+1): for m in range(N+1):
d.append([]) d.append([])
for n in range(N+1): for n in range(N+1):
d[m].append([]) d[m].append([])
for j in range(N+1): for j in range(N+1):
d[m][n].append(sp.Rational(0)) d[m][n].append(rational0)
d[0][0][0] = sp.Rational(1) d[0][0][0] = rational1
c1 = sp.Rational(1) c1 = rational1
for n in range(1, N+1): for n in range(1, N+1):
c2 = sp.Rational(1) c2 = rational1
for j in range(n): for j in range(n):
c3 = a[n] - a[j] c3 = a[n] - a[j]
c2 = c2*c3 c2 = c2*c3
for m in range(n+1):
d[m][n][j] = ((a[n] - x)*d[m][n-1][j] - m*d[m-1][n-1][j]) / c3
for m in range(n+1): for m in range(n+1):
d[m][n][n] = (c1 / c2)*(m*d[m-1][n-1][n-1] - (a[n-1] - x)*d[m][n-1][n-1]) d[m][n][j] = ((a[n] - x)*d[m][n-1][j] - m*d[m-1][n-1][j]) / c3
c1 = c2 for m in range(n+1):
coeffs = [] d[m][n][n] = (c1 / c2)*(m*d[m-1][n-1][n-1] - (a[n-1] - x)*d[m][n-1][n-1])
for m in range(len(d)): c1 = c2
coeffs.append([]) coeffs = []
for j in range(len(d)): for m in range(len(d)):
coeffs[-1].append(d[m][N][j]) coeffs.append([])
return np.array(coeffs).astype(np.float) for j in range(len(d)):
except ImportError: coeffs[-1].append(d[m][N][j])
print('didn\'t find sympy.') return np.array(coeffs).astype(np.float)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment