数值验证特定迭代法收敛性及函数选择咨询（Python/Matlab）
Absolutely, using functions of the form $f(x) = g(x)^{h(x)}$ where $g(x)=(x-r)^n$ is a fantastic choice for testing this iteration method—this structure lets you easily control root multiplicity and tweak the function's behavior to make both the math and coding straightforward. Let's walk through exactly how to pull this off, step by step.
Your proposed form is ideal because:
- You can directly set $n$ to control root multiplicity: use $n=1$ for a single root, $n\geq2$ for multiple roots (e.g., $n=2$ for a double root).
- Choosing $h(x)$ as a simple, non-zero function near $x=r$ (like a constant, $h(x)=1$, or $h(x)=x+1$) avoids introducing extra roots while letting you test both polynomial and non-polynomial cases.
- Derivatives can be computed cleanly with logarithmic differentiation, which is easy to implement symbolically (no messy manual calculus!).
We'll split this into testing single roots and multiple roots, using symbolic math to avoid derivative errors and numerical plotting to verify convergence order.
2.1 Set Up Test Functions
Single Root Example
Let’s pick $r=2$ (root location), $n=1$ (single root), and $h(x)=1$ (simple constant) for simplicity:
$$f(x) = (x-2)^1 = x-2$$
For a non-polynomial twist, use $h(x)=x+1$:
$$f(x) = (x-2)^{x+1}$$
(Note: Stick to $x>2$ for this non-polynomial case to avoid undefined logarithms.)
Multiple Root Example
Pick $r=3$ (root location), $n=2$ (double root), $h(x)=1$:
$$f(x) = (x-3)^2$$
2.2 Compute Derivatives Symbolically
Instead of calculating derivatives by hand, use sympy to generate exact expressions, then convert them to numerical functions for iteration:
import sympy as sp import numpy as np import matplotlib.pyplot as plt # Define symbolic variable x_sym = sp.symbols('x') # ---------------------- Single Root Setup ---------------------- r_single = 2 f_single = (x_sym - r_single)**1 # Simple linear function (single root) # For non-polynomial: f_single = (x_sym - r_single)**(x_sym + 1) # Compute derivatives f_prime_single = sp.diff(f_single, x_sym) f_double_prime_single = sp.diff(f_prime_single, x_sym) # Convert to numerical functions f_num_single = sp.lambdify(x_sym, f_single, 'numpy') f_p_num_single = sp.lambdify(x_sym, f_prime_single, 'numpy') f_pp_num_single = sp.lambdify(x_sym, f_double_prime_single, 'numpy') # ---------------------- Multiple Root Setup ---------------------- r_multi = 3 f_multi = (x_sym - r_multi)**2 # Double root polynomial f_prime_multi = sp.diff(f_multi, x_sym) f_double_prime_multi = sp.diff(f_prime_multi, x_sym) f_num_multi = sp.lambdify(x_sym, f_multi, 'numpy') f_p_num_multi = sp.lambdify(x_sym, f_prime_multi, 'numpy') f_pp_num_multi = sp.lambdify(x_sym, f_double_prime_multi, 'numpy')
2.3 Implement the Iteration Method
Write a function to run the iteration, track errors, and stop when convergence is reached:
def run_iteration(x0, root, f, f_p, f_pp, max_iter=50, tol=1e-12): x = x0 errors = [] for _ in range(max_iter): fx = f(x) fpx = f_p(x) fppx = f_pp(x) # Avoid division by zero/negative sqrt (start close to root!) denominator = np.sqrt(fpx**2 - fx / fppx) if denominator == 0 or np.isnan(denominator): print("Warning: Invalid denominator, stopping early.") break x_new = x - fx / denominator error = abs(x_new - root) errors.append(error) if error < tol: break x = x_new return x, errors
2.4 Verify Convergence Order
To confirm cubic convergence (single root) or linear convergence (multiple root):
- Plot error vs. iteration on a log scale (errors should drop much faster for cubic convergence).
- Calculate the slope of $\ln(e_{n+1})$ vs. $\ln(e_n)$: slopes near 3 mean cubic convergence, slopes near 1 mean linear convergence.
# Run single root iteration x0_single = 3.0 # Initial guess near root=2 root_single, errors_single = run_iteration(x0_single, r_single, f_num_single, f_p_num_single, f_pp_num_single) # Run multiple root iteration x0_multi = 4.0 # Initial guess near root=3 root_multi, errors_multi = run_iteration(x0_multi, r_multi, f_num_multi, f_p_num_multi, f_pp_num_multi) # Plot single root results plt.figure(figsize=(12, 5)) plt.subplot(121) plt.plot(errors_single, marker='o', color='b') plt.title('Single Root: Error vs Iteration') plt.xlabel('Iteration') plt.ylabel('$|x_n - r|$') plt.yscale('log') # Calculate convergence slope for single root log_errors = np.log(errors_single[:-1]) log_next_errors = np.log(errors_single[1:]) slopes_single = (log_next_errors[1:] - log_next_errors[:-1]) / (log_errors[1:] - log_errors[:-1]) plt.subplot(122) plt.plot(slopes_single, marker='s', color='b') plt.axhline(y=3, color='r', linestyle='--', label='Expected Cubic (3)') plt.title('Single Root: Convergence Order Slope') plt.xlabel('Iteration') plt.ylabel('Slope ($\ln(e_{n+1})/\ln(e_n)$)') plt.legend() # Plot multiple root results plt.figure(figsize=(12, 5)) plt.subplot(121) plt.plot(errors_multi, marker='o', color='g') plt.title('Multiple Root: Error vs Iteration') plt.xlabel('Iteration') plt.ylabel('$|x_n - r|$') plt.yscale('log') # Calculate convergence slope for multiple root log_errors_multi = np.log(errors_multi[:-1]) log_next_errors_multi = np.log(errors_multi[1:]) slopes_multi = (log_next_errors_multi[1:] - log_next_errors_multi[:-1]) / (log_errors_multi[1:] - log_errors_multi[:-1]) plt.subplot(122) plt.plot(slopes_multi, marker='s', color='g') plt.axhline(y=1, color='r', linestyle='--', label='Expected Linear (1)') plt.title('Multiple Root: Convergence Order Slope') plt.xlabel('Iteration') plt.ylabel('Slope ($\ln(e_{n+1})/\ln(e_n)$)') plt.legend() plt.show()
- Start close to the root: The iteration may diverge if your initial guess is too far away, especially for non-polynomial functions.
- Test non-polynomial cases: Swap in $h(x)=x+1$ or $h(x)=\ln(x+1)$ to confirm the method behaves consistently across function types.
- Check for numerical stability: For multiple roots, the denominator in the iteration can become sensitive—ensure your tolerance is tight enough but not so strict that you hit floating-point limits.
内容的提问来源于stack exchange，提问作者james black
