Created
August 20, 2024 19:00
-
-
Save Rodbourn/a9f4a9a1c339d997ec4dfdd7686f5a1f to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| from scipy.sparse import diags | |
| from scipy.sparse.linalg import spsolve | |
| def finite_difference_fit(x, y, order=2): | |
| """ | |
| Implement an arbitrary order one-dimensional finite difference scheme | |
| to fit a given set of data points. | |
| Parameters: | |
| x (array): x-coordinates of data points | |
| y (array): y-coordinates of data points | |
| order (int): Order of the finite difference scheme (default: 2) | |
| Returns: | |
| array: Fitted y values | |
| """ | |
| n = len(x) | |
| # Calculate step size (assuming uniform grid) | |
| h = x[1] - x[0] | |
| # Create finite difference coefficients | |
| coeffs = np.zeros(order + 1) | |
| for i in range(order + 1): | |
| coeffs[i] = (-1)**(i) * np.math.comb(order, i) | |
| # Create sparse matrix for finite difference operator | |
| diagonals = [coeffs[i] * np.ones(n - i) for i in range(order + 1)] | |
| offsets = list(range(order + 1)) | |
| A = diags(diagonals, offsets, shape=(n, n)) | |
| # Apply boundary conditions (zero derivative at boundaries) | |
| A[0, :order+1] = 0 | |
| A[0, 0] = 1 | |
| A[-1, -order-1:] = 0 | |
| A[-1, -1] = 1 | |
| # Solve the system | |
| b = y.copy() | |
| b[0] = y[0] | |
| b[-1] = y[-1] | |
| y_fit = spsolve(A, b) | |
| return y_fit | |
| # Example usage | |
| if __name__ == "__main__": | |
| # Generate some sample data | |
| x = np.linspace(0, 10, 100) | |
| y = np.sin(x) + 0.1 * np.random.randn(len(x)) | |
| # Fit using 2nd order finite difference | |
| y_fit = finite_difference_fit(x, y, order=2) | |
| # Plot results | |
| import matplotlib.pyplot as plt | |
| plt.figure(figsize=(10, 6)) | |
| plt.plot(x, y, 'o', label='Data') | |
| plt.plot(x, y_fit, label='Fitted') | |
| plt.legend() | |
| plt.xlabel('x') | |
| plt.ylabel('y') | |
| plt.title('1D Finite Difference Fit') | |
| plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment