The Stehfest approximation of the inverse transformation is a weighted mean of the solution in the Laplace domain for a number of values of the Laplace parameter $p_k$ ($k=1,...,M$). Hence, the application of the Stehfest approximation requires evaluation of $\\bar{h}(x,p)$ for a number of values of the Laplace parameter $p$ to compute the value of the head at location $x$ and time $t$. The number of terms in the approximation is $M$, which must be even. The Stehfest formula is\n", "\n", "$$\n", "h(x, t) = \\frac{\\ln 2}{t} \\sum_{k=1}^M V_k \\bar{h}(x, p_k)\n", "$$(eq_stehfest_5.4)\n", "\n", "where $M$ is the number of terms, $V_k$ is the weight for term $k$, and $p_k$ is the Laplace parameter $k$ computed as \n", "\n", "$$\n", "p_k = \\frac{k\\ln 2}{t}\n", "$$(eq_pk_5.4)\n", "\n", "In theory, the approximation is more accurate when more terms are used, but machine accuracy is quickly causing a problem. The summation in Eq. ({eq}`eq_stehfest_5.4`) consists of the summation of increasingly larger positive and negative numbers. For large values of the index $k$ in the summation, these numbers are so large that they only differ beyond the 16$^\\text{th}$ digit. Python normally only keeps track of 16 digits, which means that the sum of a positive and negative number that only differ beyond the 16$^\\text{th}$ digit essentially results in garbage. Luckily, it is fairly obvious when too many terms are used, as the result is obviously wrong.\n", "
The formula for the parameters $V_k$ is rather lengthy\n", "\n", "$$\n", "V_k = (-1)^{k+M/2} \\sum_{j=\\text{int}[(k+1)/2)]} ^ {\\min(k,M/2)}\n", "\\frac{j^{M/2}(2j)!}\n", "{(\\tfrac{M}{2}-j)!j!(j-1)!(k-j)!(2j-k)!}\n", "$$(eq_vk_5.4)\n", "\n", "The formula for $V_k$ is implemented in the routine below (adopted from the `mpmath` package by {cite:alp}`johansson2013`)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# coefficients V_k\n", "from scipy.special import factorial as fac\n", "def stehfest_coef(M):\n", " assert M % 2 == 0, 'M must be even' # make sure M is even\n", " M2 = int(M / 2)\n", " V = np.zeros(M)\n", " for k in range(1, M + 1):\n", " z = np.zeros(min(k, M2) + 1)\n", " for j in range(int((k + 1) / 2), min(k, M2) + 1):\n", " z[j] = j ** M2 * fac(2 * j) / \\\n", " (fac(M2 - j) * fac(j) * fac(j - 1) * fac(k - j) * fac(2 * j - k))\n", " V[k - 1] = (-1) ** (k + M2) * np.sum(z)\n", " return V" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ " The Stehfest algorithm ({eq}`eq_stehfest_5.4`) is implemented in the routine below. The first two input arguments are the $x$ locations and times $t$ for which the function is computed. Both input arguments `x` and `t` may be arrays. The next input argument is a Python implementation of the Laplace-transformed function that can be called as `func(p, x, keyword_arguments)`, where `keyword_arguments` can be an arbitrary number of keyword arguments. The next input argument is the keyword argument $M=12$, which is the default number of terms used in the Stehfest algorithm. Finally, an arbitrary number of keyword arguments may be supplied. These keyword arguments are passed to the function `func` of the Laplace-transformed solution (the `keyword_arguments` mentioned above). The `stehfest` function returns a two-dimensional array with shape `(len(x), len(t))`. Note that the parameters $V_k$ are only a function of $M$, so they are computed only once. The Laplace parameters $p_k$ must be computed separately for each time with Eq. ({eq}`eq_pk_5.4`)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# Stehfest algorithm\n", "def stehfest(x, t, func, M=12, **kwargs):\n", " t = np.atleast_1d(t)\n", " x = np.atleast_1d(x)\n", " f = np.zeros((len(x), len(t)))\n", " V = stehfest_coef(M)\n", " for i in range(len(t)):\n", " p = np.arange(1, M + 1) * np.log(2) / t[i]\n", " for j in range(len(x)):\n", " fbar = func(p, x[j], **kwargs)\n", " f[j, i] = np.log(2) / t[i] * np.sum(V * fbar)\n", " return f" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ " In the code blocks that follow, the transient problem of {numref}`content:chap5:step_changes_surface` (the groundwater response to a change $\\Delta h$ in the river stage in a semi-infinite confined aquifer) is re-solved using Laplace transforms. The differential equation in the Laplace domain is ({eq}`eq_laplace_transflow_5.4`), and the boundary conditions are\n", "\n", "$$\n", "\\begin{split}\n", "&\\bar{h}\\vert_{x=0}=\\frac{\\Delta h}{p} \\\\\n", "&\\bar{h}\\vert_{x\\to\\infty}=0\n", "\\end{split}\n", "$$(eq_bound_5.4)\n", "\n", "Mathematically, this is the same as the problem solved in {numref}`content:chap2:canal_to_drained_area`. The solution is \n", "\n", "$$\n", "\\bar{h} = \\frac{\\Delta h}{p} \\exp(-x\\sqrt{pS/T}) \n", "$$(eq_hbar_5.4)\n", "\n", "The solution in the time domain is obtained with the Stehfest algorithm, and the Stehfest solution is compared to the exact solution." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# parameters\n", "T = 100 # transmissivity, m^2/d\n", "S = 0.2 # storage coefficient, -\n", "delh = 2 # increase in river level, m" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# analytic solution from Section 5.1\n", "from scipy.special import erfc\n", "def h_edelman(x, t, T, S, delh=1, t0=0):\n", " u = np.sqrt(S * x ** 2 / (4 * T * (t - t0)))\n", " return delh * erfc(u)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# Laplace-transformed solution\n", "def hbar(p, x, T=100, S=1e-3, delh=2):\n", " return delh / p * np.exp(-x * np.sqrt(p * S / T))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact head at x=50 m and t=10 d: 1.2341501549039475\n", "Stehfest head with M=12 terms : [[1.23415675]]\n", "Relative error: [[5.34718718e-06]]\n" ] } ], "source": [ "# comparison\n", "x = 50 # location, m\n", "t = 10 # time, d\n", "M = 12 # number of terms\n", "V = stehfest_coef(M)\n", "hx = h_edelman(x, t, T, S, delh=2) # from Section 5.1\n", "hs = stehfest(x, t, hbar, M=M, T=T, S=S)\n", "print(f'Exact head at x={x} m and t={t} d: {hx}')\n", "print(f'Stehfest head with M={M} terms : {hs}')\n", "print(f'Relative error: {(hs - hx) / hx}')" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ " The relative error is on the order of $10^{−6}$ with $M=12$ terms in the Stehfest solution. The relative error is plotted versus the number of Stehfest terms in {numref}`Figure %s
\n",
"![](\"figures_chap5/sol_ex_5_8.png\")
\n",
"\n",
"
\n",
"\n",
"```"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.7"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": false,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Solution
\n", "\n", "$h=0.618$ m at $t=0.01$ d\\\n", "$h=1.356$ m at $t=0.1$ d\\\n", "$h=1.458$ m at $t=1$ d\n", "\n", "\n",
"\n",
"