Computing the effective resistance between nodes in finite discrete resistor grids is a classical problem in circuit analysis with applications in VLSI power delivery network analysis, graph theory, and network science. Recent advances, particularly the infinity mirror technique, provide an elegant physical interpretation for boundary conditions in finite grids. Building upon this foundation, this paper presents a closed-form analytical expression that avoids numerical truncation or polynomial fitting. Our theoretical development proceeds in two steps. First, we derive an exact analytical primitive for the singular integral term $R_2$ within the integral operator $Ω_α$. Second, we transform the doubly infinite mirror series into a compact expression using the Jacobi theta function $\vartheta_1$. This transformation achieves machine precision with only a few terms. However, under high anisotropy, the pure analytical approximation exhibits a distinct "cross-shaped" residual error. To address this, we introduce a hybrid engineering remediation: a dynamic numerical cache that performs localized grid integration (LGI), combining $O(1)$ speed with exact near-field accuracy. Numerical experiments demonstrate mean relative errors below 0.04% compared to SPICE simulations, eliminating axis-localized error artifacts. To facilitate further research, the implementation of our proposed 2D resistor grid calculator is available at: https://github.com/SeaTheDestiny/2D-Resistor-Grid-Calculator.git.