如何用Python基于收敛数值数据估计收敛结果？——针对x_n与y_n数组对的技术问询

Given your constraints (non-overlapping x arrays, increasing density, suspected power-law convergence, limited data pairs), here's a practical approach using standard Python libraries like numpy, scipy, and matplotlib:

Core Idea #

Since your data converges to a limit ( y_\infty(x) ) (likely via a power law ( y_n(x) = y_\infty(x) + C(x)n^{-\alpha} )), we can:

1. Create a target dense x-grid (for plotting) based on your largest/densest ( x_n ).
2. For each point in this grid, interpolate ( y_n ) values from all available ( (x_n, y_n) ) pairs.
3. Extrapolate these ( y_n ) sequences to ( n \to \infty ) using the suspected power-law convergence model.

Step-by-Step Implementation #

1. Setup & Sample Data (simulating your case)

First, let's simulate data matching your description (you'll replace this with your actual data):

import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Simulate converging data: y∞(x) = sin(x), y_n(x) = sin(x) + 0.5*n^(-1.2)*x + noise
np.random.seed(42)
data = []
for n in [1, 2, 3, 4]:
    x = np.linspace(0, 2*np.pi, num=10*n)  # Denser x as n increases
    y = np.sin(x) + 0.5*(n**(-1.2))*x + np.random.normal(0, 0.01, size=len(x))
    data.append({"n": n, "x": x, "y": y})

2. Create Target X-Grid

Use your densest ( x_n ) as the target grid (or create a custom linspace):

# Use the densest x array (largest n) for plotting
target_x = data[-1]["x"]

3. Estimate Global Convergence Rate (( \alpha ))

First, find the power-law exponent ( \alpha ) by analyzing how the difference between consecutive ( y_n ) decays with ( n ):

# Calculate mean absolute differences between consecutive y arrays
diffs = []
ns_diff = []
for i in range(len(data)-1):
    # Interpolate y_i to match y_{i+1}'s x grid
    y_prev_interp = interp1d(data[i]["x"], data[i]["y"], bounds_error=False, fill_value=np.nan)(data[i+1]["x"])
    mean_diff = np.nanmean(np.abs(y_prev_interp - data[i+1]["y"]))
    diffs.append(mean_diff)
    ns_diff.append(data[i]["n"])

# Fit difference ~ K*n^(-α) (linearize via log-log plot)
log_ns = np.log(ns_diff)
log_diffs = np.log(diffs)
alpha_est, log_K = np.polyfit(log_ns, log_diffs, 1)
alpha_est = -alpha_est  # Slope is -α
print(f"Estimated convergence rate α: {alpha_est:.2f}")

4. Extrapolate to Converged ( y_\infty(x) )

For each point in the target grid, extrapolate the ( y_n ) sequence to ( n \to \infty ):

def convergence_model(n, y_inf, C, alpha):
    """Power-law convergence model: y(n) = y∞ + C*n^(-α)"""
    return y_inf + C*(n**(-alpha))

target_y = []
for x_val in target_x:
    # Collect y values from all n at this x (via interpolation)
    ys = []
    ns = []
    for d in data:
        f = interp1d(d["x"], d["y"], bounds_error=False, fill_value=np.nan)
        y_val = f(x_val)
        if not np.isnan(y_val):
            ys.append(y_val)
            ns.append(d["n"])
    
    # Handle different numbers of data points
    if len(ns) == 1:
        target_y.append(ys[0])
    elif len(ns) == 2:
        # Use Richardson extrapolation with estimated α
        n1, n2 = ns
        y1, y2 = ys
        ratio = (n2/n1)**alpha_est
        y_inf = (y2*ratio - y1)/(ratio - 1)
        target_y.append(y_inf)
    else:
        # Fit the convergence model directly
        try:
            popt, _ = curve_fit(convergence_model, ns, ys, p0=[np.mean(ys), 1, alpha_est])
            target_y.append(popt[0])
        except RuntimeError:
            # Fallback to Richardson if fit fails
            n_last, n_prev = ns[-2:]
            y_last, y_prev = ys[-2:]
            ratio = (n_last/n_prev)**alpha_est
            y_inf = (y_last*ratio - y_prev)/(ratio - 1)
            target_y.append(y_inf)

target_y = np.array(target_y)

5. Plot the Result

plt.figure(figsize=(10, 6))
# Plot original data points
for d in data:
    plt.scatter(d["x"], d["y"], label=f"n={d['n']}", s=10, alpha=0.6)
# Plot converged estimate
plt.plot(target_x, target_y, "r-", linewidth=2, label="Converged Estimate")
# Plot true y∞ (for simulation validation)
plt.plot(target_x, np.sin(target_x), "k--", label="True y∞")

plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("Estimated Converged Result")
plt.show()

Key Tools & Notes #

• scipy.interpolate.interp1d: Handles interpolation between non-overlapping ( x_n ) arrays.
• scipy.optimize.curve_fit: Fits the power-law convergence model to extrapolate to ( n \to \infty ).
• Richardson Extrapolation: A robust fallback for small datasets (only 2 pairs needed if ( \alpha ) is known).

Caveats #

• If you have only 2 data pairs, the estimate depends heavily on the assumed power-law form—validate with your visualization intuition.
• Edge points (outside the overlap of all ( x_n )) may be less reliable; consider clamping your target grid to the common x-range.
• If noise is significant, use robust fitting (e.g., scipy.optimize.least_squares with a Huber loss) to reduce its impact.

内容的提问来源于stack exchange，提问作者lightfield