Poisson Jacobi

Numerical Computation Euler Methods

数学模型#

二维Poisson方程可以被表述为: $$-\nabla^2 u(x,y)=f(x,y),(x,y)\in \Omega$$

展开成坐标形式: $$-\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)=f(x,y)$$

通常使用的是Dirichlet边界条件。

网格离散#

在规则网格上使用二阶差分

$$\frac{\partial^2 u}{\partial x^2}\approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_x^2} $$

以及

$$\frac{\partial^2 u}{\partial y^2}\approx \frac{u_{i+1,j}-2u_{i,j}+u_{i+1,j}}{h_y^2} $$

带入poisson方程，有

$$-\left(\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_x^2}+\frac{u_{i+1,j}-2u_{i,j}+u_{i+1,j}}{h_y^2}\right)=f_{i,j}$$

对于$u_{i,j}$项进行收集

$$\left(\frac{2}{h_x^2}+\frac{2}{h_y^2}\right)u_{i,j}-\frac{u_{i,j+1}+u_{i,j-1}}{h_x^2}-\frac{u_{i+1,j}+u_{i-1,j}}{h_y^2}=f_{i,j}$$

可以构成一个线性系统 $Au=b$

使用Jacobi迭代求解问题#

（bug）

cpp Code#

void PoissonJacobi::solve(Field& x, const Field& b) {
    auto& u = x.data();
    const auto& rhs = b.data();

    int ny = u.rows();
    int nx = u.cols();

    Eigen::MatrixXd u_new = u;

    double hx2 = 1.0;
    double hy2 = 1.0;
    double inv = 1.0 / (2.0 / hx2 + 2.0 / hy2);

    for (int it = 0; it < max_iter_; ++it) {
        double err = 0.0;

        for (int i = 1; i < ny - 1; ++i)
            for (int j = 1; j < nx - 1; ++j) {
                double val =
                    (u(i+1,j) + u(i-1,j)) / hy2 +
                    (u(i,j+1) + u(i,j-1)) / hx2 -
                    rhs(i,j);
                u_new(i,j) = inv * val;
                err += std::abs(u_new(i,j) - u(i,j));
            }

        u.swap(u_new);
        if (err < tol_) break;
    }
}