用MATLAB对一维泊松问题 d^2u/(dx)^2=exp(-x^2)用有限差分法计算其数值解。边界条件为u(-1)=u(1)=0,-1<=x<=1

在MATLAB中，我们可以使用有限差分方法来求解一维泊松方程 \( \frac{d^2u}{dx^2} = f(x) \)，这里给定的是 \( f(x) = e^{-x^2} \)，边界条件为 \( u(-1) = u(1) = 0 \)。我们将区间 [-1, 1] 离散化成网格点，并利用中央差分法来近似导数。

以下是简单的步骤：

1. 定义网格点间隔 `h` 和计算网格点的数量 `N`（通常取 `N+1` 因为我们需要包含边界点）。
2. 创建网格点向量 `x`。
3. 对函数 \( e^{-x^2} \) 进行离散化，得到对应的函数值 `f`。
4. 初始化二维数组 `U` 来存储每个节点的数值解，边界条件初始化为0。
5. 使用有限差分公式计算中心导数 \( \frac{du}{dx} \) 和 \( \frac{d^2u}{dx^2} \)。
6. 写一个迭代循环，直到达到预设的收敛条件或者达到最大迭代次数。在循环里，更新矩阵形式下的方程，然后使用如Jacobi、Gauss-Seidel或Sor等方法解线性系统。
7. 输出最终解 `U`。

下面是一个简化的MATLAB代码示例：

function U = poisson_FD(f, N, h)
    % Define grid points and function values
    x = -1 + (0:N) * h;
    f_values = exp(-x.^2);

    % Initialize the solution array with boundary conditions
    U = zeros(N+1, 1);
    U(1) = 0;     % u(-1) = 0
    U(end) = 0;    % u(1) = 0

    % Finite difference coefficients
    D1 = -2 / h;
    D2 = 1 / (h^2);

    % Iterative loop for solving the linear system
    max_iters = 1000;
    tolerance = 1e-6;
    for iter = 1:max_iters
        Du = (D2 * U(2:end-1)) - (D1 * U(1:end-2));
        U_new = U - h^2 * Du ./ f_values;
        
        % Check convergence
        if norm(U_new - U, 'inf') < tolerance
            break;
        end
        U = U_new;
    end
end

% Call the function with your function f and desired number of nodes
f_x = @(x) exp(-x.^2);   % define the exponential function
[N, h] = deal(100);       % choose a suitable resolution
U_numerical = poisson_FD(f_x, N, h);