{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Fitting non-linear profiles\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom scipy import optimize, integrate\nimport matplotlib.pyplot as plt\nfrom multidiff import fit_crank, compute_profile_crank\nfrom derivative import dxdt\n\n\n# Define a diffusion profile from a diffusion couple,\n# where the diffusion coefficient depends exponentially on concentration\nD = .1 # diffusion coefficient in the most mobile phase\nbeta = 2.5 # coefficient of the exponential dependence\n# The diffusivity in the least mobile phase is D exp(-2 beta)\n\nl = 400\nx = np.linspace(-1, 1, l)\n\n# Concentration profile\nc = compute_profile_crank(beta, D, x)\nc += 0.03 * np.random.randn(l)\nc_interval = np.linspace(0, 1, 100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There exists a closed-form of the diffusion equation for the specific case of\n\n- a diffusion-couple geometry (two semi-infinite media of constant concentration put\n in contact at t=0)\n- an exponential dependence of the diffusivity with concentration\nTherefore we can directly fit the diffusion profile for the ``beta`` and ``D`` parameters,\nby using the function ``fit_crank``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "params = fit_crank(x, c, (1, 1))\nprint(params)\nc_fitted = compute_profile_crank(*params, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the general case of a non-constant diffusivity, the Boltzmann-Matano methods\nestimates locally the diffusivity. This methods is explained for example in \nhttps://en.wikipedia.org/wiki/Boltzmann%E2%80%93Matano_analysis \nThe local diffusion coefficient is given by\n\n\\begin{align}D(c^*) = \\frac{-1}{2t} \\frac{\\int_{c^*}^{c_L} x \\mathrm{d}c}{(\\mathrm{d}c / \\mathrm{d}x)_{x=x^*}}\\end{align}\n\nThe main difficulty is to estimate the derivative of ``c`` with ``x`` from a noisy\nsignal. It is therefore essential to smooth the signal before derivating it.\nHere we use the smooth derivative function from the ``derivative`` package. It uses\na Savitsky-Golay filter, which fits a polynomial (here of order 4) inside a local\nwindow of width 0.1. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "smooth_derivative = dxdt(c, x, kind=\"savitzky_golay\", left=.05, right=.05, order=4)\n\n# Compute the integral. In Matano's formula it is possible to compute it \n# from the left or from the right. Here there are more variations in the right\n# part of the signal so we compute the integral from the right (and hence reverse\n# the signal for the integral).\nintegrate_signal = integrate.cumtrapz(x[::-1], c[::-1])[::-1]\n\n# Matano's formula\nD_matano = -1/2 * integrate_signal / smooth_derivative[:-1]\n\nmask = np.logical_and(c > 0.25, c < 0.75)[:-1]\n\ndef exp_fit(c, bet, A):\n return A * np.exp(-2 * bet * c)\n\nparams_exp = optimize.curve_fit(exp_fit, c[:-1][mask], \n D_matano[mask])[0]\n\nprint(params_exp)\n\nfig, ax = plt.subplots(nrows=2)\nax[0].plot(x, c, 'o')\nax[0].plot(x, c_fitted)\nax[0].set_xlabel('$x$')\nax[0].set_ylabel('$C$')\nax[1].semilogy(c[:-1], D_matano, 'o', label='Matano')\nax[1].semilogy(c[:-1][mask], D_matano[mask], 'o')\nax[1].set_xlabel('$C$')\nax[1].set_ylabel('$D(C)$')\nax[1].semilogy(c_interval, D * np.exp(-2 * beta * c_interval), label='ground truth')\nax[1].semilogy(c_interval, params[1] * np.exp(-2 * params[0] * c_interval), label='fit')\nax[1].legend()\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Robustness of the fit\n\nAs we see in the figure above, the Matano estimation is noisier than the direct fit. \nIn the figure below, we compute the same estimations but for different realizations \nof the noise in the concentration profile. For some realizations, the Matano estimate\nis quite bad, while the direct fit is always very close to the ground truth.\n\nWhen a parametrization of the diffusivity is known and it is possible to compute \nnumerically the solution of the nonlinear diffusion equation with such a parametrization,\nit is therefore much robust to fit the parameters directly. However, the Matano method\nhas the great advantage of being completely generic. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=3, ncols=4, figsize=(12, 8))\nfor i in range(12):\n row = i // 4\n col = i % 4\n c = compute_profile_crank(beta, D, x)\n c += 0.03 * np.random.randn(l)\n params = fit_crank(x, c, (1, 1))\n c_fitted = compute_profile_crank(*params, x)\n\n # Matano\n smooth_derivative = dxdt(c, x, kind=\"savitzky_golay\", left=.1, right=.1, order=3)\n\n # Matano method\n integrate_signal = integrate.cumtrapz(x[::-1], c[::-1])[::-1]\n D_matano = -1/2 * integrate_signal / smooth_derivative[:-1]\n\n params_exp = optimize.curve_fit(exp_fit, c[:-1][mask], \n D_matano[mask])[0]\n print(params_exp)\n ax[row, col].semilogy(c[:-1], D_matano, 'o', label='Matano')\n ax[row, col].semilogy(c[:-1][mask], D_matano[mask], 'o')\n ax[row, col].semilogy(c_interval, D * np.exp(-2 * beta * c_interval), label='ground truth')\n ax[row, col].semilogy(c_interval, params[1] * np.exp(-2 * params[0] * c_interval), label='fit')\n ax[row, col].set_xlabel('$C$')\n ax[row, col].set_ylabel('$D(C)$')\nplt.tight_layout()\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.13" } }, "nbformat": 4, "nbformat_minor": 0 }