基于牛顿法求解Lambert W函数对应方程的收敛性分析技术问询

Let's break this down into actionable steps—first, we'll formalize the Newton-Raphson setup, handle initial guess selection, write reusable code, and generate the required convergence tables for (c=0.5, 1, 10).

1. Newton-Raphson Derivation for Our Equation #

We need to solve (f(x) = c - xe^x = 0). First, calculate the derivative of (f(x)):
[f'(x) = -(e^x + xe^x) = -e^x(1+x)]

Substitute into the standard Newton-Raphson iteration formula to get:
[x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n + \frac{c - x_n e^{x_n}}{ex(1+x_n)}]

To avoid numerical overflow for larger (x) values (like when (c=10)), we can rewrite this formula to use (e^{-x_n}) instead:
[x_{n+1} = x_n + \frac{c e^{-x_n} - x_n}{1 + x_n}]
This form stays numerically stable even as (x) grows.

2. Initial Guess Selection #

Since (xe^x) is strictly increasing for (x > -1), and our target (c) values are positive, there's exactly one positive solution for each (c>0). We can pick initial guesses based on rough estimates:

• For (c=0.5): Start with (x_0=0) (since (0*e^0=0 < 0.5))
• For (c=1): Start with (x_0=0.5) (since (0.5e^{0.5}≈0.824 < 1))
• For (c=10): Start with (x_0=2) (since (2e^2≈14.778 >10); Newton converges fast regardless of starting slightly above/below the solution)

3. Python Code Implementation #

We'll use scipy.special.lambertw to get the exact reference solution (since (W(c)) satisfies (W(c)e^{W(c)}=c)). Here's the code to run iterations and generate tables:

import math
from scipy.special import lambertw

def newton_xeq_c(c, x0, eps=1e-14, max_iter=50):
    """
    Solve xe^x = c using Newton-Raphson method.
    Returns: iteration history and exact Lambert W solution
    """
    exact_sol = lambertw(c).real
    history = []
    x_prev = x0
    
    for n in range(max_iter):
        fx = c - x_prev * math.exp(x_prev)
        abs_error = abs(x_prev - exact_sol)
        abs_fx = abs(fx)
        
        history.append((n, x_prev, abs_error, abs_fx))
        
        # Stop if both error metrics meet the threshold
        if abs_error < eps and abs_fx < eps:
            break
        
        # Compute next iteration
        exp_x = math.exp(x_prev)
        f_prime = -(exp_x + x_prev * exp_x)
        x_next = x_prev - fx / f_prime
        x_prev = x_next
    
    return history, exact_sol

# Generate tables for each c value
c_list = [0.5, 1, 10]
initial_guesses = [0, 0.5, 2]

for c, x0 in zip(c_list, initial_guesses):
    print(f"\n### Convergence Table for c={c}")
    iter_history, exact = newton_xeq_c(c, x0)
    print(f"Exact Lambert W Solution: {exact:.16f}")
    print("| Iteration | x_n               | |x_n - W(c)|       | |c - x_n e^x_n|   |")
    print("|-----------|-------------------|-------------------|-------------------|")
    for entry in iter_history:
        n, xn, err, fx_err = entry
        print(f"| {n:9d} | {xn:.16f} | {err:.16e} | {fx_err:.16e} |")

4. Convergence Tables #

Convergence Table for c=0.5 #

Exact Lambert W Solution: 0.3517337112491958

Convergence Table for c=1 #

Exact Lambert W Solution: 0.5671432904097838

Convergence Table for c=10 #

Exact Lambert W Solution: 1.7455280027406994

Key Observations #

• Newton-Raphson converges quadratically here, which is why errors plummet to (10^{-14}) in just 4-5 iterations.
• For (c > 0), we only need one solution (since (xe^x) is monotonic increasing for (x > -1)). If handling (c) in ([-1/e, 0)), we'd need two initial guesses (one below -1, one between -1 and 0) to capture both solutions.
• The rewritten iteration formula prevents overflow for large (c) values, where (e^x) could become unmanageable in the original form.

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