Python编写的数值格式中离散化误差放大问题

Hey there, let’s work through why your discretization error is blowing up for that steady-state 1D diffusion equation. I’ve tackled similar finite difference scheme issues before, so here’s a breakdown of the key checks and fixes tailored to your setup:

Your exact solution is ( u(x) = \sin(\pi x) ), so let’s plug this back into the steady-state diffusion equation ( -\frac{d^2u}{dx2} = f(x) ) (since diffusion coefficient ( D=1 )).

Calculating the second derivative: ( \frac{d^2u}{dx2} = -\pi^2 \sin(\pi x) ), which means the RHS should be ( f(x) = \pi^2 \sin(\pi x) ). If you set ( f(x) ) to something else (like just ( \sin(\pi x) )), your numerical scheme will be solving the wrong equation, leading to massive error amplification.

For ( u(x) = \sin(\pi x) ), the Dirichlet boundary conditions at ( x=0 ) and ( x=1 ) are ( u(0)=0 ) and ( u(1)=0 ). If you used Neumann (derivative) boundaries or incorrect values here, the numerical solution will diverge from the exact solution, often with growing errors.

This is a super common pitfall in Python 2.7 that breaks finite difference calculations: integer division. For example, if you calculate grid spacing as ( h = 1/N ) where ( N ) is an integer, Python 2.7 will return 0 instead of a float (since both operands are integers). This makes all your finite difference coefficients infinite or zero, which definitely causes error blowup.

Fix this by either:

• Adding from __future__ import division at the top of your script to enable float division by default
• Explicitly casting to float: ( h = float(1)/N )

For the equation ( -\frac{d^2u}{dx2} = f(x) ), the central finite difference discretization at interior points ( i ) is:
[ \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2} = f_i ]
Rearranged for the linear system ( A\mathbf{u} = \mathbf{b} ), each interior row of ( A ) should be ( [1/h^2, -2/h^2, 1/h^2] ), and ( b_i = f_i ).

If you mixed up signs (e.g., using ( -1/h^2 ) for off-diagonals) or miscalculated coefficients, the system will be ill-conditioned, leading to unstable solutions.

If you’re implementing your own iterative solver (like Jacobi or Gauss-Seidel), it might not be converging properly (e.g., wrong relaxation factor, insufficient iterations). For small to medium-sized grids, use NumPy’s built-in np.linalg.solve() instead—it’s optimized for dense linear systems and avoids convergence issues with hand-written solvers.

Example Fixed Code (Python 2.7 Compatible) #

Here’s a corrected implementation that should give minimal discretization error:

from __future__ import division  # Critical for Python 2.7 float division
import numpy as np
import matplotlib.pyplot as plt

# Grid setup
N = 100  # Number of grid points
x = np.linspace(0, 1, N)
h = x[1] - x[0]  # Safe way to calculate h (avoids integer division)

# Exact solution
u_exact = np.sin(np.pi * x)

# Correct RHS: f(x) = pi² sin(pi x)
f = (np.pi ** 2) * np.sin(np.pi * x)

# Assemble linear system A and b
A = np.zeros((N, N))
b = np.zeros(N)

# Dirichlet boundary conditions: u(0)=0, u(1)=0
A[0, 0] = 1.0
b[0] = 0.0
A[-1, -1] = 1.0
b[-1] = 0.0

# Fill interior points
for i in range(1, N-1):
    A[i, i-1] = 1 / h**2
    A[i, i] = -2 / h**2
    A[i, i+1] = 1 / h**2
    b[i] = f[i]  # Matches the discretized equation

# Solve the system
u_num = np.linalg.solve(A, b)

# Calculate and print error
max_error = np.max(np.abs(u_num - u_exact))
print("Maximum discretization error:", max_error)

# Plot results
plt.plot(x, u_exact, label="Exact Solution")
plt.plot(x, u_num, label="Numerical Solution", linestyle="--")
plt.xlabel("x")
plt.ylabel("u(x)")
plt.legend()
plt.show()

If you’re still seeing errors, try these quick checks:

• Print the value of ( h ) to confirm it’s a non-zero float
• Inspect the first few rows of matrix ( A ) to ensure coefficients are correct
• Compare numerical and exact solutions at a few key points (e.g., ( x=0.5 )) to see where discrepancies start

内容的提问来源于stack exchange，提问作者Chack.Flack