Skip to content

Instantly share code, notes, and snippets.

@kwinkunks
Created December 13, 2023 11:55
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save kwinkunks/6401dd08444da97331ad0188501be957 to your computer and use it in GitHub Desktop.
Save kwinkunks/6401dd08444da97331ad0188501be957 to your computer and use it in GitHub Desktop.
Bengt Fornbeg weights for finite difference
import numpy as np
def weights(z, x, m=0):
"""
Fornberg finite difference weights.
F90: https://github.com/bjodah/finitediff/blob/master/src/finitediff_fort.f90
Made this for Advent of Code 2023, Day 9.
Arguments:
z: location where approximations are to be accurate
x: grid point locations, found in x(0:n)
m: highest derivative for which weights are sought
Returns:
array
References:
B Fornberg (1988). Generation of Finite Difference Formulas on
Arbitrarily Spaced Grids, Mathematics of Compuation 51 (184),
p 699-706, doi: 10.1090/S0025-5718-1988-0935077-0.
Example:
>>> weights(6, range(6), 0)
array([[ -1],
[ 6],
[-15],
[ 20],
[-15],
[ 6]])
"""
n = len(x) - 1
c = np.zeros((n+1, m+1))
c1 = 1
c4 = x[0] - z
c[0, 0] = 1
for i in range(1, n+1):
mn = min(i, m)
c2 = 1
c5 = c4
c4 = x[i] - z
for j in range(i):
c3 = x[i] - x[j]
c2 *= c3
if j == i - 1:
for k in range(mn, 0, -1):
c[i, k] = c1 * (k * c[i-1, k-1] - c5 * c[i-1, k]) / c2
c[i, 0] = -c1 * c5 * c[i-1, 0] / c2
for k in range(mn, 0, -1):
c[j, k] = (c4 * c[j, k] - k * c[j, k-1]) / c3
c[j, 0] = c4 * c[j, 0] / c3
c1 = c2
return c
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment