diff --git a/book/1_gradient_divergence_curl/coord_sys.md b/book/1_gradient_divergence_curl/coord_sys.md index 4d3ccca..0caf67a 100644 --- a/book/1_gradient_divergence_curl/coord_sys.md +++ b/book/1_gradient_divergence_curl/coord_sys.md @@ -29,6 +29,18 @@ $$ \boldsymbol r = x\hat{\boldsymbol x} + y\hat{\boldsymbol y} + z\hat{\boldsymbol z}. $$ +Each base vector can be obtained by differentiating the position vector with respect to the corresponding coordinate, + +$$ +\hat{\boldsymbol x} = \partial_x\boldsymbol r, +\qquad +\hat{\boldsymbol y} = \partial_y\boldsymbol r, +\qquad +\hat{\boldsymbol z} = \partial_z\boldsymbol r. +$$ (eq:cartbase) + +We use this construction again below to find the base vectors of the cylindrical and spherical coordinate systems, where the derivatives do not come out with unit length and a normalisation factor is needed. + It is customary in the Earth Sciences to let the positive $\hat{\boldsymbol z}$-axis point downward, as depicted in {numref}`fig-cartframe`. ```{figure} figures/Cartframe.png @@ -66,7 +78,7 @@ $$ \boldsymbol r' = \hat{\boldsymbol e}_x(x-\lambda_x) + \hat{\boldsymbol e}_y(y-\lambda_y) + \hat{\boldsymbol e}_z(z-\lambda_z). $$ -Since distance should be invariant under rotation, only orthogonal tensors $\{\hat{\boldsymbol e}_x,\hat{\boldsymbol e}_y,\hat{\boldsymbol e}_z\}$ describe a rotation, and when these vectors have unit amplitude they span a base for the new frame of reference, in which case they are orthonormal vectors. +Since distance should be invariant under rotation, only orthogonal tensors $\{\hat{\boldsymbol e}_x,\hat{\boldsymbol e}_y,\hat{\boldsymbol e}_z\}$ describe a rotation, and when these vectors have unit amplitude they span a base for the new frame of reference, in which case they are orthonormal vectors, as shown in {numref}`fig-rotcartframe`. ```{figure} figures/rotCartframe.png :name: fig-rotcartframe @@ -114,6 +126,18 @@ $$ \hat{\boldsymbol x} = \hat{\boldsymbol y}\times\hat{\boldsymbol z}. $$ +You will also encounter triple vector products of the kind $\boldsymbol a\times(\boldsymbol b\times\boldsymbol c)$, which can be evaluated as + +$$ +\boldsymbol a\times(\boldsymbol b\times\boldsymbol c) = (\boldsymbol a\cdot\boldsymbol c)\boldsymbol b - (\boldsymbol a\cdot\boldsymbol b)\boldsymbol c . +$$ (eq:tripprod) + +The result lies in the plane spanned by $\boldsymbol b$ and $\boldsymbol c$, because $\boldsymbol b\times\boldsymbol c$ is perpendicular to that plane and the second cross product turns the result back into it. The order of the brackets matters, $\boldsymbol a\times(\boldsymbol b\times\boldsymbol c)\ne(\boldsymbol a\times\boldsymbol b)\times\boldsymbol c$. The following identity can be useful as well, + +$$ +\boldsymbol a\times(\boldsymbol b\times\boldsymbol c) + \boldsymbol b\times(\boldsymbol c\times\boldsymbol a) = (\boldsymbol a\times\boldsymbol b)\times\boldsymbol c . +$$ (eq:tripprod2) + ## Curvilinear coordinate systems Apart from the Cartesian reference frame, which is a rectangular coordinate system, two coordinate systems are often employed that use curved base vectors. One is the cylindrical coordinate system and the other is the spherical coordinate system. The point shown in {numref}`fig-cartframe` is depicted in cylindrical and spherical coordinate systems in {numref}`fig-cylsphere`. @@ -282,13 +306,26 @@ $$ ## Exercises 1. What is the length of the vector function $\boldsymbol v$? -2. Show that $\hat{\boldsymbol x} = \partial\boldsymbol r/\partial x$. -3. A two-dimensional rectangle in the plane $z=0$ is spanned by the two vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$ and $\boldsymbol d_y = d_y\hat{\boldsymbol y}$, where $d_x$ and $d_y$ are the lengths in the $x$- and $y$-directions, respectively, and $\hat{\boldsymbol n}$ is the unit normal on the rectangle pointing in the positive $z$-direction. Give a geometrical interpretation of the product $\hat{\boldsymbol n}\cdot(\boldsymbol d_x\times\boldsymbol d_y)$. -4. A three-dimensional rectangle, a brick, with dimensions $d_x,d_y,d_z$ is spanned by the three vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$, $\boldsymbol d_y = d_y\hat{\boldsymbol y}$ and $\boldsymbol d_z = d_z\hat{\boldsymbol z}$. Give a geometrical interpretation of the product $\boldsymbol d_x\cdot(\boldsymbol d_y\times\boldsymbol d_z)$. -5. Consider a smooth surface $\mathbb{S}$ with unique unit normal vector $\hat{\boldsymbol n}$. Show that any vector quantity $\boldsymbol H$ can be composed as $\boldsymbol H = (\hat{\boldsymbol n}\cdot\boldsymbol H)\hat{\boldsymbol n} + (\hat{\boldsymbol n}\times\boldsymbol H)\times\hat{\boldsymbol n}$ and give a geometric interpretation of the two terms. -6. A smooth surface $\mathbb{S}$ has a unique unit normal vector $\hat{\boldsymbol n}$. Decompose the gradient operator into a part tangential and a part normal to the surface $\mathbb{S}$. -7. Use the same reasoning as in the exercise above and find $\hat{\boldsymbol r} = A_r\,\partial\boldsymbol r/\partial r$, $\hat{\boldsymbol\phi} = A_\phi\,\partial\boldsymbol r/\partial\phi$, and $\hat{\boldsymbol\theta} = A_\theta\,\partial\boldsymbol r/\partial\theta$. Carry out the differentiations and find the coefficients $(A_r,A_\phi,A_\theta)$ by requiring that these vectors have unit length. These results should lead to an expression for the unit vectors of the spherical coordinate system in Cartesian coordinates. -8. The upper right-hand side of the matrix of {numref}`tab-direction-cosines` is left empty. The $3\times3$ block in the bottom-left corner of the table relates the unit vectors of the spherical reference frame to those of the Cartesian reference frame. Suppose we write this as +2. Show that $\boldsymbol a\cdot(\boldsymbol b\times\boldsymbol c) = \boldsymbol b\cdot(\boldsymbol c\times\boldsymbol a) = \boldsymbol c\cdot(\boldsymbol a\times\boldsymbol b)$. +3. Evaluate + + $$ + \begin{aligned} + \hat{\boldsymbol x}&\cdot(\hat{\boldsymbol y}\times\hat{\boldsymbol z}), &\qquad + \hat{\boldsymbol x}&\cdot(\hat{\boldsymbol z}\times\hat{\boldsymbol y}), \\ + \hat{\boldsymbol y}&\cdot(\hat{\boldsymbol x}\times\hat{\boldsymbol z}), &\qquad + \hat{\boldsymbol y}&\cdot(\hat{\boldsymbol z}\times\hat{\boldsymbol x}), \\ + \hat{\boldsymbol z}&\cdot(\hat{\boldsymbol y}\times\hat{\boldsymbol x}), &\qquad + \hat{\boldsymbol z}&\cdot(\hat{\boldsymbol x}\times\hat{\boldsymbol y}). + \end{aligned} + $$ +4. A two-dimensional rectangle in the plane $z=0$ is spanned by the two vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$ and $\boldsymbol d_y = d_y\hat{\boldsymbol y}$, where $d_x$ and $d_y$ are the lengths in the $x$- and $y$-directions, respectively, and $\hat{\boldsymbol n}$ is the unit normal on the rectangle pointing in the positive $z$-direction. Give a geometrical interpretation of the product $\hat{\boldsymbol n}\cdot(\boldsymbol d_x\times\boldsymbol d_y)$. +5. A three-dimensional rectangle, a brick, with dimensions $d_x,d_y,d_z$ is spanned by the three vectors $\boldsymbol d_x = d_x\hat{\boldsymbol x}$, $\boldsymbol d_y = d_y\hat{\boldsymbol y}$ and $\boldsymbol d_z = d_z\hat{\boldsymbol z}$. Give a geometrical interpretation of the product $\boldsymbol d_x\cdot(\boldsymbol d_y\times\boldsymbol d_z)$. +6. Evaluate $(\boldsymbol a\times\boldsymbol b)\times\boldsymbol c$. +7. Consider a smooth surface $\mathbb{S}$ with unique unit normal vector $\hat{\boldsymbol n}$. Show that any vector quantity $\boldsymbol H$ can be composed as $\boldsymbol H = (\hat{\boldsymbol n}\cdot\boldsymbol H)\hat{\boldsymbol n} + (\hat{\boldsymbol n}\times\boldsymbol H)\times\hat{\boldsymbol n}$ and give a geometric interpretation of the two terms. +8. A smooth surface $\mathbb{S}$ has a unique unit normal vector $\hat{\boldsymbol n}$. Decompose the gradient operator into a part tangential and a part normal to the surface $\mathbb{S}$. +9. Find $\hat{\boldsymbol r} = A_r\,\partial\boldsymbol r/\partial r$, $\hat{\boldsymbol\phi} = A_\phi\,\partial\boldsymbol r/\partial\phi$, and $\hat{\boldsymbol\theta} = A_\theta\,\partial\boldsymbol r/\partial\theta$. Carry out the differentiations and find the coefficients $(A_r,A_\phi,A_\theta)$ by requiring that these vectors have unit length. These results should lead to an expression for the unit vectors of the spherical coordinate system in Cartesian coordinates. +10. The upper right-hand side of the matrix of {numref}`tab-direction-cosines` is left empty. The $3\times3$ block in the bottom-left corner of the table relates the unit vectors of the spherical reference frame to those of the Cartesian reference frame. Suppose we write this as $$ \left(\begin{array}{c} \hat{\boldsymbol r} \\ \hat{\boldsymbol\phi} \\ \hat{\boldsymbol\theta}\end{array}\right) @@ -301,7 +338,7 @@ $$ -\sin(\phi) & \cos(\phi) & 0 \\ \cos(\phi)\cos(\theta) & \sin(\phi)\cos(\theta) & -\sin(\theta) \end{array}\right). - $$ + $$ (eq:Rmat) Show that @@ -311,5 +348,16 @@ $$ $$ where $\mathsf{R}^{t}$ denotes the transpose of $\mathsf{R}$. Then show that $\mathsf{R}\mathsf{R}^{t} = \mathsf{I}$, where $\mathsf{I}$ is the $3\times3$ unit matrix, implying that the transpose of $\mathsf{R}$ is equal to its inverse. This follows from the fact that $\mathsf{R}$ is a matrix of orthonormal vectors and is therefore an orthonormal matrix. -9. Give explicit expressions for the base vectors of the spherical reference frame by carrying out the matrix-vector multiplication of {eq}`eq:sphCar`. Notice that all three base vectors depend on both angles. This implies that the base vectors are not constant but depend on position! The spherical reference frame is not an inertial frame. -10. Carry out the following differentiations: $\partial\hat{\boldsymbol r}/\partial\theta$, $\partial\hat{\boldsymbol\theta}/\partial\theta$, $\partial\hat{\boldsymbol\phi}/\partial\theta$ and $\partial\hat{\boldsymbol r}/\partial\phi$, $\partial\hat{\boldsymbol\theta}/\partial\phi$, $\partial\hat{\boldsymbol\phi}/\partial\phi$. +11. Give explicit expressions for the base vectors of the spherical reference frame by carrying out the matrix-vector multiplication of {eq}`eq:sphCar`, with the rotation matrix given by {eq}`eq:Rmat`. You should find that + + $$ + \hat{\boldsymbol r} = \left(\begin{array}{c} \cos(\phi)\sin(\theta) \\ \sin(\phi)\sin(\theta) \\ \cos(\theta)\end{array}\right), + \qquad + \hat{\boldsymbol\theta} = \left(\begin{array}{c} \cos(\phi)\cos(\theta) \\ \sin(\phi)\cos(\theta) \\ -\sin(\theta)\end{array}\right), + \qquad + \hat{\boldsymbol\phi} = \left(\begin{array}{c} -\sin(\phi) \\ \cos(\phi) \\ 0\end{array}\right). + $$ + + Notice that each base vector depends on at least one angle. This implies that the base vectors are not constant but depend on position! The spherical reference frame is not an inertial frame. +12. Evaluate $\hat{\boldsymbol r}\times\hat{\boldsymbol\phi}$, $\hat{\boldsymbol r}\times\hat{\boldsymbol\theta}$, and $\hat{\boldsymbol\phi}\times\hat{\boldsymbol\theta}$. +13. Carry out the following differentiations: $\partial\hat{\boldsymbol r}/\partial\theta$, $\partial\hat{\boldsymbol\theta}/\partial\theta$, $\partial\hat{\boldsymbol\phi}/\partial\theta$ and $\partial\hat{\boldsymbol r}/\partial\phi$, $\partial\hat{\boldsymbol\theta}/\partial\phi$, $\partial\hat{\boldsymbol\phi}/\partial\phi$, and express them all in terms of $\hat{\boldsymbol r}$, $\hat{\boldsymbol\theta}$ and $\hat{\boldsymbol\phi}$. diff --git a/book/1_gradient_divergence_curl/curl.md b/book/1_gradient_divergence_curl/curl.md index d73a939..e763c67 100644 --- a/book/1_gradient_divergence_curl/curl.md +++ b/book/1_gradient_divergence_curl/curl.md @@ -82,22 +82,112 @@ $$ \oint_{\boldsymbol r}\boldsymbol\tau\cdot\boldsymbol v\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}\hat{\boldsymbol n}\cdot(\nabla\times\boldsymbol v)\,\mathrm{d}S . $$ (eq:stokes) +## Vector identities with the curl + +The curl as expressed in {eq}`eq:curl` can be written as a matrix-vector operation, + +$$ +\nabla\times\boldsymbol v = \left(\begin{array}{ccc} +0 & -\partial_z & \partial_y \\ +\partial_z & 0 & -\partial_x \\ +-\partial_y & \partial_x & 0 +\end{array}\right)\left(\begin{array}{c} v_x \\ v_y \\ v_z\end{array}\right). +$$ (eq:curlmat) + +The matrix is antisymmetric and its diagonal is empty, which is another way of saying that no component of $\boldsymbol v$ contributes to the component of $\nabla\times\boldsymbol v$ along its own direction. + +The curl of a curl is given by + +$$ +\nabla\times\nabla\times\boldsymbol v = \nabla(\nabla\cdot\boldsymbol v) - \nabla^2\boldsymbol v, +$$ (eq:curlcurl) + +where $\nabla^2 = \nabla\cdot\nabla$ is the Laplacian. The curl of the curl of a vector field is the gradient of the divergence of that field, given by the first term on the right-hand side, minus the divergence of the gradient of that field, given by the second term. In this case the evaluation in terms of vector operations is the same as in {eq}`eq:tripprod`, with $\boldsymbol a$ and $\boldsymbol b$ both replaced by $\nabla$. + +If the curl is taken of a vector product of two vector functions, the product rule for differentiation must be taken into account as well, + +$$ +\nabla\times(\boldsymbol v\times\boldsymbol w) = \boldsymbol v(\nabla\cdot\boldsymbol w) - \boldsymbol w(\nabla\cdot\boldsymbol v) + (\boldsymbol w\cdot\nabla)\boldsymbol v - (\boldsymbol v\cdot\nabla)\boldsymbol w . +$$ (eq:curlvecprod) + +The first two terms take the divergence of $\boldsymbol w$ and of $\boldsymbol v$ and point the result along $\boldsymbol v$ and $-\boldsymbol w$ respectively, while the last two terms take the derivative of $\boldsymbol v$ and of $\boldsymbol w$ along the direction of $\boldsymbol w$ and $-\boldsymbol v$ respectively. + +The divergence operator can also act on the vector product of two functions, $\nabla\cdot(\boldsymbol v\times\boldsymbol w)$. This is evaluated with the vector rules of the exercises on coordinate systems, keeping the product rule for differentiation in mind, because the $\nabla$-operator acts on both $\boldsymbol v$ and $\boldsymbol w$. It is given by + +$$ +\nabla\cdot(\boldsymbol v\times\boldsymbol w) = \boldsymbol w\cdot(\nabla\times\boldsymbol v) - \boldsymbol v\cdot(\nabla\times\boldsymbol w). +$$ (eq:divvecprod) + +The divergence of a vector product of two functions takes the curl of each of them, computes the scalar product with the other, and subtracts the results. + +## The magnetic field around an infinite current carrying wire + +In the chapter on Potential Fields you will investigate this problem a bit deeper, but here we use it to illustrate that a field that circulates around another field does not have a non-zero curl everywhere, and that it can have zero divergence as well. Let us take an infinitely long current carrying wire along the vertical axis. That makes it a two-dimensional problem, because nothing is changing in the $z$-direction. In cylindrical coordinates the magnetic field is given by + +$$ +\boldsymbol H = \frac{I}{2\pi\varrho}\hat{\boldsymbol\phi}, +$$ (eq:HJcyl) + +which in Cartesian coordinates is written as + +$$ +\boldsymbol H = \frac{I}{2\pi\varrho^2}\left(-y\,\hat{\boldsymbol x} + x\,\hat{\boldsymbol y}\right). +$$ (eq:HJcart) + +The curl of this magnetic field has only a vertical component, as you can imagine from Ampère's law, because it states that $\nabla\times\boldsymbol H = \boldsymbol J$, with $\boldsymbol J$ the electric current density in the wire. We have to look only at the vertical component, because $\boldsymbol J$ has only a vertical component. We find + +$$ +\hat{\boldsymbol z}\cdot(\nabla\times\boldsymbol H) = \partial_x H_y - \partial_y H_x = \frac{I}{2\pi}\left(\partial_y\frac{y}{\varrho^2} + \partial_x\frac{x}{\varrho^2}\right) = 0, \qquad \text{for } \varrho>0 . +$$ + +We find that the magnetic field is curl-free everywhere outside the wire. The field rotates around the wire but it has zero curl. If we take the divergence of the magnetic field, we find + +$$ +\nabla\cdot\boldsymbol H = \partial_x H_x + \partial_y H_y = \frac{I}{2\pi}\left(-\partial_x\frac{y}{\varrho^2} + \partial_y\frac{x}{\varrho^2}\right) = 0, \qquad \text{for } \varrho>0 . +$$ + +We find that the magnetic field is divergence free as well, everywhere outside the wire. We investigate what happens at the wire. For that we take a planar circular area in the $(x,y)$-plane, as a cross-section that goes through the wire, but we do not specify the radius. We multiply Ampère's law with a unit vector in the vertical direction and integrate over the surface area $\mathbb{S}$ to find + +$$ +\int_{\mathbb{S}}\hat{\boldsymbol z}\cdot(\nabla\times\boldsymbol H)\,\mathrm{d}x\,\mathrm{d}y = \int_{\mathbb{S}}\hat{\boldsymbol z}\cdot\boldsymbol J\,\mathrm{d}x\,\mathrm{d}y . +$$ + +The right-hand side is identified as the current $I$ running in the wire, and we use Stokes' theorem on the left-hand side to arrive at + +$$ +\oint_{\partial\mathbb{S}}\boldsymbol\tau\cdot\boldsymbol H\,\mathrm{d}l . +$$ + +This is a line integral around the full circle, which is the boundary of the circular surface area that we started out with, and $\boldsymbol\tau$ is the unit tangent vector. In cylindrical coordinates, $\boldsymbol\tau\,\mathrm{d}l = \hat{\boldsymbol\phi}\,\varrho\,\mathrm{d}\phi$. Substituting this together with the expression of the magnetic field in {eq}`eq:HJcyl` in the above equation results in + +$$ +\oint_{\partial\mathbb{S}}\boldsymbol\tau\cdot\boldsymbol H\,\mathrm{d}l = \frac{I}{2\pi}\int_{\phi=0}^{2\pi}\mathrm{d}\phi = I, +$$ + +which is equal to the surface integral with the normal component of the current. In conclusion, the magnetic field has zero curl at any point outside the wire, but the line integral around any circle around the wire of the tangential component of the magnetic field is equal to the current flowing in the wire. This will come back when we discuss Potential Fields, and again when we discuss Electromagnetic Fields. + +## Sources of non-zero curl + +Rigid rotation and shear are two known causes for generating a field with non-zero curl. + ## Exercises 1. Use Stokes' theorem to evaluate $\oint_{\boldsymbol r}\boldsymbol\tau\,\mathrm{d}l$, where $\boldsymbol\tau$ is the unit tangent along the closed boundary of the area $\mathbb{S}$. The integration runs in the direction of circulation that forms a right-handed system with the unit normal vector on $\mathbb{S}$. -2. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\,p(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $p(\boldsymbol r)$ is a continuously differentiable scalar function, Stokes' integral theorem gives +2. Use the matrix expression for the curl in {eq}`eq:curlmat` to find the matrix expression for the curl of the curl of $\boldsymbol v$, as expressed in vector form in {eq}`eq:curlcurl`. +3. Give a detailed derivation to show that {eq}`eq:divvecprod` is correct. +4. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\,p(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $p(\boldsymbol r)$ is a continuously differentiable scalar function, Stokes' integral theorem gives $$ \oint_{\boldsymbol r}\boldsymbol\tau\,p\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}(\hat{\boldsymbol n}\times\nabla)p\,\mathrm{d}S, $$ which is Stokes' theorem for the gradient. -3. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\times\boldsymbol w(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $\boldsymbol w(\boldsymbol r)$ is a continuously differentiable vector function, Stokes' integral theorem gives +5. Show that when $\boldsymbol v(\boldsymbol r) = \boldsymbol a\times\boldsymbol w(\boldsymbol r)$, where $\boldsymbol a$ is an arbitrary constant vector and $\boldsymbol w(\boldsymbol r)$ is a continuously differentiable vector function, Stokes' integral theorem gives $$ \oint_{\boldsymbol r}\boldsymbol\tau\times\boldsymbol w\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}(\hat{\boldsymbol n}\times\nabla)\times\boldsymbol w\,\mathrm{d}S . $$ -4. Ampère's law is given by $\nabla\times\boldsymbol H = \boldsymbol J$, which states that the electric current is equal to the curl of the magnetic field. Convert this to integral form using Stokes' theorem, by using a flat surface $\mathbb{S}$ and choosing a unit normal vector $\hat{\boldsymbol n}$ on $\mathbb{S}$. You should find +6. Ampère's law is given by $\nabla\times\boldsymbol H = \boldsymbol J$, which states that the electric current is equal to the curl of the magnetic field. Convert this to integral form using Stokes' theorem, by using a flat surface $\mathbb{S}$ and choosing a unit normal vector $\hat{\boldsymbol n}$ on $\mathbb{S}$. You should find $$ \oint_{\boldsymbol r}\boldsymbol H\cdot\boldsymbol\tau\,\mathrm{d}l = \int_{\boldsymbol r\in\mathbb{S}}\hat{\boldsymbol n}\cdot\boldsymbol J\,\mathrm{d}S . diff --git a/book/1_gradient_divergence_curl/divergence.md b/book/1_gradient_divergence_curl/divergence.md index 1e97bf2..035de39 100644 --- a/book/1_gradient_divergence_curl/divergence.md +++ b/book/1_gradient_divergence_curl/divergence.md @@ -55,6 +55,140 @@ We have just derived Gauss' integral theorem, or divergence theorem, and we foun The divergence of a vector field is the outward flux of the vector field per unit volume. ::: +### Gauss' theorem in one, two and three dimensions + +In one dimension the vector function can have only one component, say $v_x$, and it can depend on the spatial variable $x$ only. A finite volume in one dimension is a line segment, and its surface consists of the two end points of that segment. Let the line run from $x=x_b$ to $x=x_e$, from beginning to end. The divergence is $\nabla\cdot\boldsymbol v = \partial_x v_x$ and Gauss' theorem reads + +$$ +\int_{x_b}^{x_e}\partial_x v_x\,\mathrm{d}x = v_x(x_e) - v_x(x_b). +$$ (eq:gauss1d) + +You will recognise this as the fundamental theorem of calculus, stating that integration reverses differentiation. How do we know that it is the term $v_x(x_b)$ that should be given the minus sign? + +In two dimensions a vector field can have two components and depend on two spatial variables. Take the $(x,y)$-plane as our two-dimensional space. The vector field is written as $\boldsymbol v(x,y) = v_x(x,y)\hat{\boldsymbol x} + v_y(x,y)\hat{\boldsymbol y}$, and we evaluate the outward flow out of the rectangular domain $x_b0$ that is valid until the ball hits the ground again for the first time. + + (b) Determine the time $T_0$ at which the ball hits the ground and the maximum height $H$ the ball reaches the first time. + + Assume that every time the ball hits the ground there is some plastic deformation of the ball and the ball loses a fraction $\gamma$ of its energy at every bounce, such that at the $n^{\text{th}}$ bounce the velocity can be expressed as $v_n = \sqrt{1-\gamma}\,v_{n-1}$ for $n>0$. + + (c) The time it takes to complete the $n^{\text{th}}$ bounce, $T_n$, can be written in terms of $T_{n-1}$ as $T_n = \alpha T_{n-1}$, and in terms of $T_0$ as $T_n = \beta T_0$. Determine $\alpha$ and $\beta$. + + (d) Show that $T_n = (1-\gamma)^{n/2}\sqrt{\dfrac{8H}{g}}$. + + (e) Give an expression for $T_\infty$ by evaluating the series $T_\infty = \sum_{m=0}^{\infty}T_m$. Since $\gamma$ is a fraction, it has a positive non-zero value smaller than 1, so that the series can be summed to an explicit expression. Why does it take only a finite time for the ball to make infinitely many bounces? + + (f) If the energy loss per bounce is extremely small, $\gamma\ll 1$, use a Taylor series expansion around $\gamma=0$ for $T_\infty$ and keep only the leading term as a first-order approximation to the time required to make an infinite number of bounces. + + (g) Plot the height as a function of time and try different values for $\gamma$, to understand truncation of a series as an approximation of a function. How many terms do you need in the sum to arrive at the same time as the approximate expression? + + :::{admonition} Answers to (a)–(f) + :class: tip dropdown + (a) $z = v_0 t - \tfrac{1}{2}gt^2$. + + (b) $H = \dfrac{v_0^2}{2g}$ and $T_0 = \dfrac{2v_0}{g}$. + + (c) $T_n$ is proportional to $v_n$, hence $T_n = \sqrt{1-\gamma}\,T_{n-1} = (1-\gamma)^{n/2}T_0$, so $\alpha=\sqrt{1-\gamma}$ and $\beta=(1-\gamma)^{n/2}$. + + (d) In the expression for $T_n$ in (c), express $T_0$ in terms of $H$ instead of in terms of $v_0$. + + (e) $T_\infty = T_0\sum_{m=0}^{\infty}(1-\gamma)^{m/2} = \sqrt{\dfrac{8H}{g}}\dfrac{1}{1-\sqrt{1-\gamma}}$. + + (f) In that case $\sqrt{1-\gamma}\approx 1-\tfrac{\gamma}{2}$ and $T_\infty \approx \sqrt{\dfrac{8H}{g}}\dfrac{2}{\gamma}$. + ::: diff --git a/book/3_diffusion_fields/em_diffusion.md b/book/3_diffusion_fields/em_diffusion.md new file mode 100644 index 0000000..73133a8 --- /dev/null +++ b/book/3_diffusion_fields/em_diffusion.md @@ -0,0 +1,228 @@ +# Electromagnetic diffusive field equations with a source + +Maxwell's equations are given by + +$$ +-\nabla\times\boldsymbol H + \varepsilon\,\partial_t\boldsymbol E + \sigma\boldsymbol E = -\boldsymbol J^e, +$$ (eq:hrel) + +$$ +\nabla\times\boldsymbol E + \mu\,\partial_t\boldsymbol H = \boldsymbol 0, +$$ (eq:erel) + +where the right-hand sides contain the sources that generate the electromagnetic field. In {eq}`eq:hrel`, the source is the electric current dipole $\boldsymbol J^e$ (A/m$^2$), and in {eq}`eq:erel` the source is zero, because no magnetic monopoles have been found. Some electric sources can be approximated by a magnetic current dipole source. A closed current-carrying loop antenna is such an electric current source that is often described as a magnetic dipole, which can be understood from integrating both sides of {eq}`eq:erel` over the area of the loop, performing a scalar product of each term with the unit normal to the loop, and using Stokes' theorem on the curl of the electric field. In the left-hand sides, $\boldsymbol E$ and $\boldsymbol H$ are the electric and magnetic field strengths in V/m and A/m, respectively, and the medium parameters $\varepsilon,\sigma,\mu$ are the electric permittivity and conductivity, and the magnetic permeability, with units of s/($\Omega$m), 1/($\Omega$m), and s$\Omega$/m, respectively. These equations describe in principle the electric, magnetic and electromagnetic fields for all possible time variations, including the static fields discussed before. Here we are interested in the diffusive field regime. We see that {eq}`eq:hrel`, Ampère's law modified by Maxwell, contains both the electric field and the time derivative of the electric field. If we take the divergence of this equation, we obtain + +$$ +(\sigma + \varepsilon\,\partial_t)\,\nabla\cdot\boldsymbol E = -\nabla\cdot\boldsymbol J^e . +$$ + +We use the fact that $\boldsymbol D=\varepsilon\boldsymbol E$ and that $\nabla\cdot\boldsymbol D=\rho_f$, with $\rho_f$ being the volume density of free charge (C m$^{-3}$). Substituting this in the equation gives + +$$ +\partial_t\rho_f = -\nabla\cdot(\sigma\boldsymbol E + \boldsymbol J^e), +$$ + +which is known as the equation of continuity of electric current, or the charge conservation law. We replace the source by assuming there is an accumulation of free charge $\rho_0$ at $t=0$, which removes the source and replaces it with the initial condition $\rho_f(0)=\rho_0$. If we write the electric field as $\boldsymbol E=\boldsymbol D/\varepsilon$, the equation becomes + +$$ +\partial_t\rho_f = -\frac{\sigma}{\varepsilon}\rho_f, +$$ + +the solution of which is known as + +$$ +\rho_f = \rho_0\exp\left(-\frac{\sigma}{\varepsilon}t\right) = \rho_0\exp(-t/\tau_r), +$$ + +where $\tau_r=\varepsilon/\sigma$ is known as the charge relaxation time. It is a measure of the time it takes the medium to return to its equilibrium state after it has been disturbed by an electromagnetic wave. First, we learn from this equation the following theorem: *within a region of non-vanishing conductivity there can be no permanent distribution of free charge*. This is the significant equation, because when alternating source currents have oscillation periods much larger than the relaxation time, the effect of $\varepsilon$ can be neglected and we can effectively assume $\varepsilon=0$ in the first Maxwell equation. For most Earth materials, $\varepsilon\lesssim 10^{-9}$ F/m and $\sigma\gtrsim 10^{-4}$ S/m, so for most earth materials under investigation the relaxation time is less than 10 $\mu$s. Hence, for source time function oscillations with periods longer than 0.3 ms, we can safely take $\varepsilon=0$. For now, we assume this is the case. Electromagnetic waves are treated later in the course. We write Maxwell's equations in the diffusive approximation out in components as + +$$ +\begin{aligned} +-\partial_y H_z + \partial_z H_y + \sigma E_x &= -J_x^e, \\ +-\partial_z H_x + \partial_x H_z + \sigma E_y &= -J_y^e, \\ +-\partial_x H_y + \partial_y H_x + \sigma E_z &= -J_z^e, \\ +\partial_y E_z - \partial_z E_y + \mu\,\partial_t H_x &= 0, \\ +\partial_z E_x - \partial_x E_z + \mu\,\partial_t H_y &= 0, \\ +\partial_x E_y - \partial_y E_x + \mu\,\partial_t H_z &= 0 . +\end{aligned} +$$ + +First we link these equations to the 1D heat equation by assuming a one-dimensional problem here as well. If the source is an infinite current sheet with a constant current in the $(x,y)$-plane, it means that all fields are independent of the $x$- and $y$-coordinates and all derivatives $\partial_x$ and $\partial_y$ become zero. Then, when the current has an $x$-component only, we have $J_y^e=J_z^e=0$. We substitute these findings in the equations and find + +$$ +\begin{aligned} +\partial_z H_y + \sigma E_x &= -J_x^e, \\ +-\partial_z H_x + \sigma E_y &= 0, \\ +\sigma E_z &= 0, \\ +-\partial_z E_y + \mu\,\partial_t H_x &= 0, \\ +\partial_z E_x + \mu\,\partial_t H_y &= 0, \\ +\mu\,\partial_t H_z &= 0, +\end{aligned} +$$ (eq:plzHy) + +which eliminates the vertical electric and magnetic field components and results in two equal but independent sets of two coupled differential equations, one for $E_x,H_y$ and one for $E_y,H_x$, but the latter has no source to generate the fields. The source is taken as $J_x^e(z,t)=\delta(z)W(t)$, which represents the plane electric current sheet in the $(x,y)$-plane at $z=0$ with a current directed along the $x$-axis and with a time function $W(t)$. This means it is a point source in the $z$-direction and of infinite extent in the $x$- and $y$-directions. In the course Signals and Time Series, the delta-function was introduced in time, but we can use it in all coordinates. Here, this leads to equations for $E_x$ and $H_y$, as + +$$ +\begin{aligned} +\partial_z H_y(z,t) + \sigma E_x(z,t) &= -J_x^e(z,t), \\ +\partial_z E_x(z,t) + \mu\,\partial_t H_y(z,t) &= 0 . +\end{aligned} +$$ + +The diffusive field that satisfies these equations is called the Transverse ElectroMagnetic (TEM) diffusive field. The diffusion equation for the electric field is obtained by eliminating the magnetic field from the two equations above. We do that as follows, + +$$ +\begin{aligned} +\partial_z\left(\mu\,\partial_t H_y(z,t)\right) + \sigma\mu\,\partial_t E_x(z,t) &= -\mu\,\partial_t J_x^e(z,t), \\ +\mu\,\partial_t H_y(z,t) &= -\partial_z E_x(z,t), +\end{aligned} +$$ + +and substituting the right-hand side of the second equation to replace the term with the magnetic field in the first equation. The electric field is found to satisfy the following diffusion equation: + +$$ +\partial_z\partial_z E_x - \mu\sigma\,\partial_t E_x = \mu\,\partial_t J_x^e(z,t). +$$ (eq:dee) + +For $J_x^e(z,t)=\delta(z)W(t)$, the solution is known and given by + +$$ +E_x(z,t) = -\int_{t'=0}^{t}\left(\frac{\sigma\mu z^2}{4(t-t')} - \frac{1}{2}\right)\sqrt{\frac{\mu}{4\pi\sigma(t-t')^3}}\exp\left(-\frac{\sigma\mu z^2}{4(t-t')}\right)W(t')\,\mathrm{d}t', +$$ (eq:extW) + +where the convolutional integral is present because the function $W(t)$ has not been defined yet. If the current is switched on and kept at a constant value $I$ for all times after $t=0$, $W(t)=0$ for $t<0$ and $W(t)=I$ for $t>0$, the integral can be evaluated explicitly and the solution is + +$$ +E_x(z,t) = -I\sqrt{\frac{\mu}{4\pi\sigma t}}\exp\left(-\frac{\sigma\mu z^2}{4t}\right), \quad\text{for } t>0 . +$$ (eq:ExtH) + +You recognise the similarity with the solution to the heat equation with a Gaussian distribution for the initial temperature on an infinitely long wire. Note that here we make explicit that the expression for the electric field exists only for positive times. This is because $t=0$ is the moment the source is switched on, and physics tells us that the Earth will respond to a source action only at or after the moment the source becomes active, and not before. This is the notion of causality. + +## Solution using Fourier transformation + +In the course on Signals and Time Series, you have learned Fourier series and the Fourier transformation. In the heat equation for a finite-length wire, we use Fourier series to find solutions. We then showed that for an infinitely long wire that has a Gaussian temperature distribution at zero time, a solution in terms of Fourier series is possible as well. We also showed that if we let the distance between two frequency values go to zero, the Fourier transformation is found. But we can solve partial differential equations with a non-zero source term more formally with the aid of Fourier transformations. We do this by transforming {eq}`eq:dee` to the frequency domain using + +$$ +\hat{\boldsymbol E}(z,\omega) = \int_{t=0}^{\infty}\exp(-\mathrm{i}\omega t)\,\boldsymbol E(z,t)\,\mathrm{d}t, +$$ + +to obtain + +$$ +\partial_z\partial_z\hat E_x - \mathrm{i}\omega\sigma\mu\,\hat E_x = \mathrm{i}\omega\mu\,\hat W(\omega)\,\delta(z). +$$ (eq:dees) + +The diacritical hat on the quantities denotes the Fourier transform of the quantity. The time derivative is gone and replaced by an algebraic multiplication with $\mathrm{i}\omega$, and the equation has been reduced to an ordinary second-order differential equation. As you have learned in the course Signals and Time Series, the exponential function is an eigenfunction, so we can propose a solution in the form + +$$ +\hat E_x = \hat A\exp\left(-\sqrt{\mathrm{i}\omega\mu\sigma}\,|z|\right), +$$ + +where $\hat A$ is the constant that needs to be determined. The reason why $|z|$ is in the argument of the exponential function is that the field must diffuse away from the source and decrease in amplitude with increasing distance in both the positive and negative $z$-directions. Evaluating the first and second derivatives with respect to $z$ gives + +$$ +\begin{aligned} +\partial_z\hat E_x &= -\sqrt{\mathrm{i}\omega\sigma\mu}\,\hat A\exp\left(-\sqrt{\mathrm{i}\omega\sigma\mu}\,|z|\right)\frac{z}{|z|}, \\ +\partial_z\partial_z\hat E_x &= \mathrm{i}\omega\mu\sigma\,\hat A\exp\left(-\sqrt{\mathrm{i}\omega\sigma\mu}\,|z|\right) - 2\sqrt{\mathrm{i}\omega\sigma\mu}\,\hat A\,\delta(z)\exp\left(-\sqrt{\mathrm{i}\omega\sigma\mu}\,|z|\right), \\ +&= \mathrm{i}\omega\sigma\mu\,\hat E_x - 2\sqrt{\mathrm{i}\omega\sigma\mu}\,\hat A\,\delta(z), +\end{aligned} +$$ + +where we have used the fact that + +$$ +\mathrm{sign}(z) = \frac{z}{|z|} = \left\{\begin{array}{rl} -1, & z<0, \\ 1, & z>0, \end{array}\right. +$$ + +is known as the sign-function. It has a sudden step of size 2 at $z=0$, for which reason its derivative is equal to $\partial_z\,\mathrm{sign}(z)=2\delta(z)$. Because the impulse function acts only at $z=0$, we can substitute that already in the exponential function that multiplies the delta-function. If we substitute these results in {eq}`eq:dees` we find that + +$$ +\hat A = -\frac{\sqrt{\mathrm{i}\omega\mu}}{2\sqrt{\sigma}}\hat W . +$$ + +The solution for the electric field is then found as + +$$ +\hat E_x(z,\omega) = -\frac{\sqrt{\mu}}{2\sqrt{\sigma}}\hat W\sqrt{\mathrm{i}\omega}\exp\left(-\sqrt{\mathrm{i}\omega\sigma\mu}\,|z|\right). +$$ (eq:eesd) + +This is the space-frequency domain solution for the one-dimensional electric field in the Earth in response to an electric current sheet that is an impulse at $z=0$ but exists for all $x$- and $y$-values, and has a time function $W(t)$. The transformation back to time results in the expression given in {eq}`eq:extW`, and in {eq}`eq:ExtH` in case the source is a step-function switch-on, in which case its frequency transform is given by $\hat W=I/(\mathrm{i}\omega)$. Now we can explain the convolution in time, because you know from the course Signals and Time Series that a product in the frequency domain leads to a convolution in the time domain. + +## The concept of the principle of superposition + +Before we go to the three-dimensional problem, we take a look at the right-hand side of {eq}`eq:dee` for the given source time function that is switched on at $t=0$ and kept constant after that, and note that $\partial_t u(t)=\delta(t)$, where $u(t)$ is the unit-step function and $\delta(t)$ is the delta-function. In that case the equation for the electric field becomes + +$$ +\partial_z\partial_z E_x - \mu\sigma\,\partial_t E_x = \mu I\,\delta(t)\,\delta(z). +$$ (eq:deeg) + +The function that collapses a differential equation to an impulse in space and time is called the Green's function. Physically, the Green's function describes the (scalar) response to an impulse in space and time. This is a general principle for linear partial differential equations and you have discussed this in the course Signals and Time Series. The equation for the space-time impulse response is then + +$$ +(\partial_z\partial_z - \mu\sigma\,\partial_t)\,G(z-z',t-t') = -\delta(t-t')\,\delta(z-z'), +$$ (eq:deg) + +where the minus sign in the right-hand side is a historic choice. The constant displacements $t'$ and $z'$ are chosen because the real source may have a finite length in the $z$-direction and have a certain time duration. If we now multiply both sides of this equation with a function $-\mu\,\partial_{t'}J_x^e(z',t')$, assume that $J_x^e(z',t')\ne 0$ only for $z_bt', +$$ (eq:G3D) + +and it satisfies + +$$ +(\nabla^2 - \sigma\mu\,\partial_t)\,G = -\delta(\boldsymbol r-\boldsymbol r')\,\delta(t-t') . +$$ + +Again we recognise the similarity with the solution for the heat equation for an infinitely long wire with an initial temperature that was a Gaussian distribution. It also resembles the 1D diffusive field solution to Maxwell's equations. The major difference is that here the late-time behaviour is proportional to $t^{-3/2}$. The electric and magnetic fields have much more complicated behaviour, but that is far beyond the scope of the present course. + +### Brief analysis of the outcome + +For a source in the location $\boldsymbol r'$, the field diffuses away in spherical direction and the amplitude only depends on radial distance. The product $\sigma\mu$ is the diffusion constant and has the unit of s/m$^2$. It is a measure of how fast the width of the Gaussian spreads in space for a given time instant. + +**Early time behaviour.** The early-time behaviour is similar to what we saw in the 1D heat equation. + +**Late time behaviour.** The late-time behaviour is inversely proportional to $t^{3/2}$, which is faster than what we saw in the solutions to the 1D heat and Maxwell diffusion equations. + +**Time to the peak.** If we differentiate the Green's function to $t$ and equate the result to zero, we find the time to the peak to be $t-t'=t_p$, at $t_p=\sigma\mu|\boldsymbol r-\boldsymbol r'|^2/6$. The amplitude of the Green's function at the peak time is given by + +$$ +G(\boldsymbol r-\boldsymbol r',t-t'=t_p) = \frac{3\sqrt{6}\exp(-3/2)}{4\pi^{3/2}\sigma\mu|\boldsymbol r-\boldsymbol r'|^3}, +$$ + +which shows that if you double the distance to the source, the peak amplitude is 8 times smaller. diff --git a/book/3_diffusion_fields/figures/GHeat.png b/book/3_diffusion_fields/figures/GHeat.png new file mode 100644 index 0000000..b648c93 Binary files /dev/null and b/book/3_diffusion_fields/figures/GHeat.png differ diff --git a/book/3_diffusion_fields/figures/GHeatA.png b/book/3_diffusion_fields/figures/GHeatA.png new file mode 100644 index 0000000..65cf0d6 Binary files /dev/null and b/book/3_diffusion_fields/figures/GHeatA.png differ diff --git a/book/3_diffusion_fields/figures/heat_ic_quarters.svg b/book/3_diffusion_fields/figures/heat_ic_quarters.svg new file mode 100644 index 0000000..a800540 --- /dev/null +++ b/book/3_diffusion_fields/figures/heat_ic_quarters.svg @@ -0,0 +1,426 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/book/3_diffusion_fields/figures/heat_ic_thirds.svg b/book/3_diffusion_fields/figures/heat_ic_thirds.svg new file mode 100644 index 0000000..c4742c2 --- /dev/null +++ b/book/3_diffusion_fields/figures/heat_ic_thirds.svg @@ -0,0 +1,388 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/book/3_diffusion_fields/figures/make_figures.py b/book/3_diffusion_fields/figures/make_figures.py new file mode 100644 index 0000000..6f21438 --- /dev/null +++ b/book/3_diffusion_fields/figures/make_figures.py @@ -0,0 +1,43 @@ +"""Regenerate the initial-condition sketches of the heat-equation exercises. + +Run from this directory: python make_figures.py +""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +INK = "#6b7280" # neutral grey, readable on light and dark backgrounds +LINE = "#00A6D6" # TU Delft cyan + +def step_sketch(fname, edges, values, xticks, xticklabels): + fig, ax = plt.subplots(figsize=(4.2, 2.4)) + xs, ys = [], [] + for (a, b), v in zip(zip(edges[:-1], edges[1:]), values): + xs += [a, b] + ys += [v, v] + ax.plot(xs, ys, color=LINE, lw=2) + ax.axhline(0, color=INK, lw=0.8) + ax.set_xlim(-0.05, 1.12) + ax.set_ylim(-1.35, 1.35) + ax.set_xticks(xticks, xticklabels) + ax.set_yticks([-1, 0, 1]) + ax.set_xlabel("$x$", color=INK) + ax.set_ylabel("$h(x)$", color=INK) + for side in ("top", "right"): + ax.spines[side].set_visible(False) + for side in ("left", "bottom"): + ax.spines[side].set_color(INK) + ax.tick_params(colors=INK) + fig.tight_layout() + fig.savefig(fname, transparent=True, metadata={"Date": None}) + plt.close(fig) + +# exercise 4: h = 0, -1, 1, 0 on the four quarters of the bar, T(0,t) = T(L,t) = 0 +step_sketch("heat_ic_quarters.svg", + [0, 0.25, 0.5, 0.75, 1.0], [0, -1, 1, 0], + [0, 0.25, 0.5, 0.75, 1.0], ["0", "$L/4$", "$L/2$", "$3L/4$", "$L$"]) + +# exercise 5: h = 1, -1, 1 on the three thirds of the bar, insulated ends +step_sketch("heat_ic_thirds.svg", + [0, 1/3, 2/3, 1.0], [1, -1, 1], + [0, 1/3, 2/3, 1.0], ["0", "$L/3$", "$2L/3$", "$L$"]) diff --git a/book/3_diffusion_fields/heat_equation_1d.md b/book/3_diffusion_fields/heat_equation_1d.md new file mode 100644 index 0000000..4f32414 --- /dev/null +++ b/book/3_diffusion_fields/heat_equation_1d.md @@ -0,0 +1,533 @@ +# The heat equation in one dimension + +Let us assume the heat can flow only in the $x$-direction, like in a bar. This reduces the heat equation, {eq}`eq:Heat`, to + +$$ +\partial_t T - \kappa\,\partial_x^2 T = 0, +$$ (eq:Heat1D) + +where we assume no source or sink is present, for which reason the right-hand side is zero. If we want to find a solution to this equation, we will need other information. These are initial conditions and boundary conditions. The number of conditions needed is equal to the highest order derivative in a variable. {eq}`eq:Heat1D` has a first-order time derivative, which means we need one initial or boundary condition, and it has a second-order spatial derivative, which means we need two initial or boundary conditions in space. Let us assume the bar has length $L$ and extends from $x=0$ to $x=L$. So-called Dirichlet conditions describe the temperature at the end points, $x=0$ and $x=L$, e.g., $T(0,t)=0$ and $T(L,t)=0$. So-called Neumann conditions prescribe the normal component of the current, or, in this case, the first spatial derivative of temperature. For example, a no-flow end-point condition for an insulated bar, meaning no heat can flow away from the bar, is expressed as + +$$ +\begin{aligned} +\partial_x T(0,t) &= 0, \\ +\partial_x T(L,t) &= 0. +\end{aligned} +$$ + +For some source-free problems, solutions can be found using the technique of separation of variables. This means that we can write the temperature in the bar as a product of a function that depends only on position and a function that depends only on time. We are going to find a solution for two different cases. + +## Example 1: bar with zero temperature at the end points + +For such a bar we use the 1D heat equation with conditions on the temperature at the end points, + +$$ +\begin{aligned} +\partial_t T &= \kappa\,\partial_x^2 T, \\ +T(x,0) &= h(x), \\ +T(0,t) &= 0, \\ +T(L,t) &= 0. +\end{aligned} +$$ + +Note that the boundary condition at zero time is an arbitrary function $h(x)$ of position. This condition allows for non-trivial solutions and this means that the function $h(x)$ can be seen as the starting point of the diffusive process. We could see it as equivalent to the presence of a source. The boundary conditions do not depend on time and space simultaneously and we can try a solution with separation of variables. We take a solution in the form + +$$ +T(x,t) = f(x)\,g(t), +$$ + +and put it in the equation to find + +$$ +\begin{aligned} +\partial_t\left[f(x)g(t)\right] &= \kappa\,\partial_x^2\left[f(x)g(t)\right], \\ +f(x)\,\partial_t g(t) &= \kappa\,g(t)\,\partial_x^2 f(x), \\ +\frac{\partial_t g(t)}{g(t)} &= \kappa\frac{\partial_x^2 f(x)}{f(x)} = -\kappa\lambda . +\end{aligned} +$$ + +The separation constant $\lambda$ is introduced to be able to separate the two equations and solve them independently. Note that attaching $\kappa$ to the separation constant is not required but convenient. If you don't do it, it does not matter. When you discover in a later step that it would be better to do it, you can come back here and do it after all. We can now write the equation as two equations, + +$$ +\begin{aligned} +\partial_t g(t) &= -\kappa\lambda\,g(t), \\ +\partial_x^2 f(x) &= -\lambda\,f(x). +\end{aligned} +$$ + +Now, we must find out if we can satisfy the boundary conditions for $T(x,t)$ and move them to conditions for $f(x)$ and $g(t)$. We find + +$$ +T(0,t) = f(0)g(t) = 0 \quad\text{and}\quad T(L,t) = f(L)g(t) = 0 . +$$ + +If we would choose $g(t)=0$ then we would find the trivial solution that the temperature is zero everywhere and for all times. Hence, we choose $f(0)=f(L)=0$ as boundary conditions for $f(x)$. The separation constant $\lambda$ is known as the eigenvalue of the problem. We are going to find that non-trivial solutions are found only for specific values of $\lambda$. The non-trivial solution for a value of $\lambda$ is then called an eigenfunction. You have encountered or will encounter this in the course ECTB2130, Signals and Time Series. An eigenfunction of a differential equation is a function that, when it is taken as an input signal, gives as output the same function scaled by a factor, which is then the corresponding eigenvalue. + +Now we are ready to solve the corresponding two equations, + +$$ +\begin{aligned} +(\partial_t + \kappa\lambda)\,g &= 0, \\ +(\partial_x^2 + \lambda)\,f &= 0, \\ +f(0) = f(L) &= 0 . +\end{aligned} +$$ + +In fact, we solve first the spatial equation. We must consider three ranges for $\lambda$ to find solutions for the second-order differential equation. + +**Case $\lambda>0$.** For positive eigenvalues, the solution consists of oscillating functions, + +$$ +f(x) = c_1\sin(x\sqrt{\lambda}) + c_2\cos(x\sqrt{\lambda}). +$$ (eq:fsincosd) + +Because $f(0)=0$, we find $c_2=0$. Because $f(L)=0$, we find that $\sin(L\sqrt{\lambda})=0$, which leads to + +$$ +\lambda = \left(\frac{n\pi}{L}\right)^2, \quad\text{for } n=1,2,3,\cdots, +$$ + +hence the solution is given by + +$$ +f(x) = c_1\sin\left(\frac{n\pi x}{L}\right), \quad\text{for } n=1,2,3,\cdots . +$$ + +There is an infinite number of eigenvalues, each of which has its own eigenfunction as solution to the problem. Note that $c_1$ is not determined by this condition and we could remove it, because later we will need to make sure to satisfy the last condition of $T(x,0)=h(x)$. For simplicity we keep it here as an undetermined constant. + +**Case $\lambda=0$.** For a zero eigenvalue, the solution is given by $f(x)=c_0+c_1x$. Because both $f(0)$ and $f(L)$ must be zero, we find that both coefficients are zero and only the trivial solution exists. Hence $f(x)=0$. This also means that $\lambda=0$ is not an eigenvalue for this problem! + +**Case $\lambda<0$.** For negative eigenvalues, the solution consists of exponential functions given by + +$$ +f(x) = c_1\exp(x\sqrt{-\lambda}) + c_2\exp(-x\sqrt{-\lambda}). +$$ + +Setting $f(0)=0$ and $f(L)=0$ leads again to $c_1=0$ and $c_2=0$. We see that for negative values of $\lambda$ only the trivial solution is found. Hence, also no negative eigenvalues exist for this problem. + +Now it is time to solve the time equation. We must consider only positive values of $\lambda$ and substitute $\lambda_n$ that we can later specify as found above. For each value of $\lambda_n$, we should find a solution $g_n(t)$. The corresponding equation is written as + +$$ +(\partial_t + \kappa\lambda_n)\,g_n = 0 . +$$ + +This is a simple first-order equation and we know that exponential functions are repeated when differentiated. Hence, we find the solution to be of the form + +$$ +g_n(t) = c_n\exp(-\kappa\lambda_n t) = c_n\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right]. +$$ + +The temperature is given by the product of the two solutions. We can write for every $n$, + +$$ +T_n(x,t) = C_n\sin\left(\frac{n\pi x}{L}\right)\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right]. +$$ (eq:T1Dn) + +The coefficient $C_n$ combines the coefficients we found in the solution in space and in time. Note that this is a solution, but that does not mean that this function alone can satisfy all boundary conditions we can think of. The function $h(x)$ of the boundary condition at zero time is an arbitrary function of position. This means that the solution is quite restrictive and can be used only for functions $h(x)$ that match $T_n(x,0)$. We are of course allowed to combine solutions for different values of $n$, because the problem is linear and the principle of superposition applies. This means we can write the temperature as a solution to the heat equation as + +$$ +T(x,t) = \sum_{n=1}^{\infty} C_n\sin\left(\frac{n\pi x}{L}\right)\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], +$$ + +with initial condition + +$$ +T(x,0) = h(x) = \sum_{n=1}^{\infty} C_n\sin\left(\frac{n\pi x}{L}\right). +$$ (eq:BCheat1D) + +With an infinite sum of functions, we can imagine we can satisfy almost all functions $h(x)$. Now the question is: How can we find the coefficients $C_n$ such that {eq}`eq:BCheat1D` is satisfied? That is not very difficult. We see that the argument of the sine-function is at least spanning half a period for $00$.** For positive eigenvalues, the solution consists of oscillating functions, + +$$ +f(x) = c_1\sin(x\sqrt{\lambda}) + c_2\cos(x\sqrt{\lambda}). +$$ (eq:fsincosn) + +Because $f'(0)=0$, we find $\left.\partial_x f(x)\right|_{x=0} = c_1\sqrt{\lambda} = 0$, hence $c_1=0$. Because $f'(L)=0$, we find that $c_2\sqrt{\lambda}\sin(L\sqrt{\lambda})=0$, which leads to + +$$ +\lambda = \left(\frac{n\pi}{L}\right)^2, \quad\text{for } n=1,2,3,\cdots, +$$ + +hence the solution is given by + +$$ +f_n(x) = c_2\cos\left(\frac{n\pi x}{L}\right), \quad\text{for } n=1,2,\cdots . +$$ + +There is an infinite number of eigenvalues, each of which has its own eigenfunction as solution to the problem. + +**Case $\lambda=0$.** For a zero eigenvalue, the solution is given by $f(x)=c_0+c_1x$. Because both $f'(0)$ and $f'(L)$ must be zero, we find that the coefficient $c_1=0$ and only the solution $f(x)=c_0$ exists for $\lambda=0$. This means that $\lambda=0$ is an eigenvalue for this problem! + +**Case $\lambda<0$.** For negative eigenvalues, the solution consists of exponential functions given by + +$$ +f(x) = c_1\exp(x\sqrt{-\lambda}) + c_2\exp(-x\sqrt{-\lambda}). +$$ + +To impose $f'(0)=0$ and $f'(L)=0$ we need + +$$ +f'(x) = \sqrt{-\lambda}\left[c_1\exp(x\sqrt{-\lambda}) - c_2\exp(-x\sqrt{-\lambda})\right]. +$$ + +This leads from $f'(0)=0$ to $c_1=c_2$, and from $f'(L)=0$ we find $c_1=c_2=0$. Hence, there are no negative eigenvalues for this problem. + +The total solution so far is + +$$ +f_n(x) = c_2\cos\left(\frac{n\pi x}{L}\right), \quad\text{for } n=0,1,2,\cdots . +$$ + +Now it is time to solve the time equation. We recognise that this part of the problem is the same as before and we find the solution to be of the form + +$$ +g_n(t) = c_n\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right]. +$$ + +The temperature is given by the product of the two solutions. We can write for every $n$ + +$$ +T_n(x,t) = C_n\cos\left(\frac{n\pi x}{L}\right)\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], \quad n=0,1,2,\cdots . +$$ (eq:T1Dnc) + +Note that $n=0$ is now part of the solution. The temperature as a solution to the heat equation is the sum of all solutions $T_n$ and given by + +$$ +T(x,t) = \sum_{n=0}^{\infty} C_n\cos\left(\frac{n\pi x}{L}\right)\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], +$$ + +with initial condition + +$$ +T(x,0) = h(x) = \sum_{n=0}^{\infty} C_n\cos\left(\frac{n\pi x}{L}\right). +$$ + +What we found for Fourier sine-series in the previous example applies to the Fourier cosine-series as well. Applying Fourier orthogonality of cosine-functions results in the final solution + +$$ +\begin{aligned} +T(x,t) &= \sum_{n=0}^{\infty} C_n\cos\left(\frac{n\pi x}{L}\right)\exp\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], \\ +C_n &= \left\{\begin{array}{ll} +\dfrac{1}{L}\displaystyle\int_{x=0}^{L} h(x)\,\mathrm{d}x, & n=0, \\[3mm] +\dfrac{2}{L}\displaystyle\int_{x=0}^{L} h(x)\cos\left(\frac{n\pi x}{L}\right)\mathrm{d}x, & n=1,2,3,\cdots . +\end{array}\right. +\end{aligned} +$$ + +With these two expressions, we have fully solved a partial differential equation with Neumann boundary conditions. Note that now the temperature does not necessarily go to zero when time goes to infinity. + +## A completely different solution procedure + +### An infinitely long wire with an initial Gaussian temperature distribution + +Solving problems with separation of variables is not always the best solution, e.g., when the temperature distribution at $t=0$ is a Gaussian function. For this example, we assume an infinitely long wire. Hence, we are solving + +$$ +\begin{aligned} +\partial_t T &= \kappa\,\partial_x^2 T, \\ +T(x,0) &= T_0\exp\left(-\frac{x^2}{L^2}\right). +\end{aligned} +$$ (eq:c0time) + +Note that we do not use explicit boundary conditions at the end points, but are satisfied with zero conditions, which seems to make sense from a physical point of view. We can see from the initial condition of {eq}`eq:c0time` that + +$$ +\lim_{x\rightarrow\pm\infty} T(x,t) = 0 . +$$ (eq:endcond) + +The boundary condition at zero time shows that at $x=0$ and $t=0$ the temperature is $T_0$, which is the peak value of the temperature. The constant $L$ in the Gaussian is a measure of the width of the Gaussian function. The actual peak value may change and the width of the Gaussian may change, but we assume that the solution will have a Gaussian structure in space, and we are going to find what its time behaviour then is. Let us propose a solution of the form + +$$ +T(x,t) = f(t)\exp\left(-h(t)x^2\right). +$$ + +The two functions $f(t)$ and $h(t)$ allow for time variations in peak temperature and the time at which the peak occurs. Note that this solution cannot be easily written as a product of a time-dependent function and a space-dependent function as we did in the method of separation of variables. For $h(t)\ge 0$, the temperature will exponentially decrease with increasing $x$. Because we introduce two new time functions that we want to find, we need boundary conditions for them. We can find them from the boundary condition on $T(x,0)$ as + +$$ +f(0) = T_0, \qquad h(0) = \frac{1}{L^2} . +$$ (eq:initchf) + +We substitute the proposed solution in {eq}`eq:Heat1D` and find + +$$ +\begin{aligned} +\exp\left(-h(t)x^2\right)\left[\partial_t f(t) - x^2 f(t)\,\partial_t h(t)\right] &= \kappa\left[4f(t)h^2(t)x^2 - 2f(t)h(t)\right]\exp\left(-h(t)x^2\right), \\ +\partial_t f(t) + 2\kappa f(t)h(t) &= x^2\left[f(t)\,\partial_t h(t) + 4\kappa f(t)h^2(t)\right]. +\end{aligned} +$$ + +Similar to the solution with separation of variables, here we see that the left-hand side of the equation constitutes an equation that can be solved independent from $x$, and the right-hand side does depend on the location. Hence the whole equation can be satisfied for all values of $x$ only when the left-hand side is satisfied independently from the right-hand side and vice versa. We find the following two differential equations, + +$$ +\partial_t f(t) = -2\kappa f(t)h(t), +$$ (eq:eqhf) + +$$ +\partial_t h(t) = -4\kappa h^2(t). +$$ (eq:eqhh) + +The first equation depends on both $f(t)$ and $h(t)$ and must be solved after we solve the equation for $h(t)$. It is not hard to show that + +$$ +h(t) = \frac{1}{4\kappa t + L^2} +$$ + +is a solution to {eq}`eq:eqhh` that satisfies the initial condition of {eq}`eq:initchf` as well. Similarly, it is then not hard to show that + +$$ +f(t) = T_0\frac{L}{\sqrt{4\kappa t + L^2}} +$$ + +satisfies {eq}`eq:eqhf` and satisfies the initial condition of {eq}`eq:initchf`. This solves the problem and we find the temperature along the bar as a function of time as + +$$ +T(x,t) = \frac{T_0 L}{\sqrt{4\kappa t + L^2}}\exp\left(-\frac{x^2}{4\kappa t + L^2}\right). +$$ (eq:TGauss) + +Notice that we have solved the heat equation using the initial condition for $T(x,0)$ at zero time expressed in {eq}`eq:c0time` and end-point conditions at infinity for all times expressed in {eq}`eq:endcond`. We found a solution in terms of a single expression using ordinary functions. It is clear that this problem cannot be easily solved with the method of separation of variables. In the exercises, we find a solution in terms of the Fourier cosine series. The explicit exact solution that we found here gives a better insight in the space-time behaviour of diffusive fields than the solution with separation of variables, and can be analysed. + +### Analysis of the outcome + +**Early time behaviour.** We can see that for zero time, the solution is well behaved, because the square root in time in the denominator goes to $L$ when time goes to zero and the temperature goes to $T_0$ as given as the initial value. If we think about the origin of this heat, it must have come from some negative time, and we can imagine it must have come from time equal to $t=-L^2/(4\kappa)$. In that case the denominator goes to zero and the polynomial term will blow up to infinity. But we can see that for $x\ne 0$, the exponential function has the time in the denominator of its argument and that will bring the temperature to zero at $t=-L^2/(4\kappa)$. The only problematic point is then $x=0$, which must be the location of the origin of the heat. This is characteristic for solutions to the diffusive field equation with a source that is approximated to be localised in a single point. Note that all finite order derivatives to time of this function go to zero at zero time. Another way of saying this is: an exponential function with a negative argument that is inversely proportional to time goes to zero faster than a polynomial of any finite inverse power of time goes to infinity for zero time. Notice that for times $t<-L^2/(4\kappa)$ the solution is no longer physical. The reason is that the square root in the polynomial term will have a negative argument and the outcome of such a square root will be complex, which we regard as non-physical in space-time. You can see that + +$$ +\begin{aligned} +\lim_{t\uparrow -\frac{L^2}{4\kappa}}\exp\left(-\frac{x^2}{4\kappa t + L^2}\right) &= \infty, \\ +\lim_{t\rightarrow -\infty}\exp\left(-\frac{x^2}{4\kappa t + L^2}\right) &= 1, +\end{aligned} +$$ + +from which we conclude that negative infinite time is not a problem for the exponential term, but the negative side of the zero denominator in the exponential function is the problem. In addition, we can see that for any fixed time $t<-L^2/(4\kappa)$ the exponential function is no longer well-behaved as a function of $x$ and grows exponentially for increasing distance. That is a second reason to regard solutions for $t<-L^2/(4\kappa)$ as non-physical. + +**Late time behaviour.** The late-time behaviour is entirely determined by the polynomial function, and hence we conclude that the temperature drops proportional to the inverse square root of time and does not depend on location anymore! The reason is that when $t\rightarrow\infty$, $\exp\left[-x^2/(4\kappa t + L^2)\right]\rightarrow 1$ for any fixed position. Hence, the exponential function is a Gaussian in space, but as a function of time it is zero at zero time and it smoothly goes to 1 at infinite time. Hence, for late times, the exponential function has no more influence on the value of the temperature in the bar. + +**Time to the peak temperature.** We can compute how much time it takes for the temperature to reach its peak value at any point on the bar. To determine the time to the peak value, we take the derivative to time and equate the result to zero. The derivative to time of {eq}`eq:TGauss` is given by + +$$ +\begin{aligned} +\partial_t T(x,t) &= \exp\left(-\frac{x^2}{4\kappa t + L^2}\right)\partial_t\frac{T_0 L}{\sqrt{4\kappa t + L^2}} ++ \frac{T_0 L}{\sqrt{4\kappa t + L^2}}\,\partial_t\exp\left(-\frac{x^2}{4\kappa t + L^2}\right), \\ +&= -\frac{2\kappa T_0 L}{(4\kappa t + L^2)^{3/2}}\exp\left(-\frac{x^2}{4\kappa t + L^2}\right) ++ \frac{4\kappa T_0 L x^2}{(4\kappa t + L^2)^{5/2}}\exp\left(-\frac{x^2}{4\kappa t + L^2}\right). +\end{aligned} +$$ + +This derivative is zero when time is equal to the time to the peak, $t_p$, and this is given by + +$$ +t_p = \frac{2x^2 - L^2}{4\kappa}, +$$ + +from which it is clear that only for $|x|>L/\sqrt{2}$ this happens at positive times. That means that for points on the bar where $|x|L/\sqrt{2}$ the temperature will first rise for $0L/\sqrt{2} . +$$ + +We conclude that the temperature has a maximum value that is inversely proportional to distance $|x|$, and it occurs on the bar for positions $|x|>L/\sqrt{2}$. + +To show the early- and late-time behaviour, we show the temperature as a function of time. In the example, $T_0=1$, $\kappa=1$, $L=10$, and we compute the temperature for positions $-30L/\sqrt{2}\approx 7$ m occurs at later times with increasing distance to the mid point. The time of the maximum temperature value grows with distance squared. This is not very well visible in the plot, because of the logarithmic time scale. + +```{figure} figures/GHeat.png +:name: fig-gheat +:width: 75% + +The temperature as a function of space and time with a Gaussian distribution in space at $t=0$ as initial condition. +``` + +Looking at the temperature as a function of time for specific locations away from the mid point is useful to evaluate the interplay between the late-time, or polynomial, part and the early-time, or exponential, part. To this end, {numref}`fig-gheata` shows them as a function of time for distances of $x=10$ m (left plot) and $x=20$ m (right plot) as black curves, together with the polynomial part (red dashed lines) and the exponential part (blue dashed lines). The vertical thin black lines mark the time to the maximum temperature, which is 25 s for $x=10$ m and 175 s for $x=20$ m. + +```{figure} figures/GHeatA.png +:name: fig-gheata +:width: 100% + +The temperature as a function of time at $x=10$ m (left) and $x=20$ m (right), with a Gaussian distribution in space at $t=0$ as initial condition. The black curve is the temperature, the red dashed line the polynomial part, the blue dashed line the exponential part, and the thin vertical line marks the time $t_p$ to the maximum temperature. +``` + +### An infinitely long wire with an initial double Gaussian temperature distribution + +If the initial condition specifies two shifted Gaussians, the solution can be written immediately by understanding that the superposition principle holds. This means that two separate initial conditions each lead to a solution that can be summed and that is equal to the solution we would obtain by specifying the two initial conditions as a single starting point and then solving the equation. Hence, if we have + +$$ +\begin{aligned} +\partial_t T &= \kappa\,\partial_x^2 T, \\ +T(x,0) &= T_0\left[\exp\left(-\frac{(x+2L)^2}{L^2}\right) + \exp\left(-\frac{(x-2L)^2}{L^2}\right)\right], +\end{aligned} +$$ (eq:2c0time) + +the solution can be readily given by + +$$ +T(x,t) = \frac{T_0 L}{\sqrt{4\kappa t + L^2}}\left[\exp\left(-\frac{(x+2L)^2}{4\kappa t + L^2}\right) + \exp\left(-\frac{(x-2L)^2}{4\kappa t + L^2}\right)\right]. +$$ (eq:T2Gauss) + +From this solution, and the fact that the exponential functions are quite distinct for early times but become equal for late times, we can understand that after some time we cannot distinguish between a field that consisted initially of two Gaussians and one that consisted of a single Gaussian temperature distribution. The closer the initial peak values are together, the earlier the blurring is such that the distinction cannot be made anymore. In a nutshell, that is the difficulty with analysing diffusive field signals when we record them to obtain information about the Earth. + +## The infinitely long wire revisited: a solution with the Fourier cosine series + +Solving problems with separation of variables is not always the best solution, for example when the temperature distribution at $t=0$ is a Gaussian function, but it is possible and can still be convenient. For the example, we return to the infinitely long wire. Hence, we are solving the same equation and boundary condition as before, + +$$ +\begin{aligned} +\partial_t T &= \kappa\,\partial_x^2 T, \\ +T(x,0) &= T_0\exp\left(-\frac{x^2}{L^2}\right). +\end{aligned} +$$ + +This function goes to zero at the end points "infinity", which seems to make sense from a physical point of view. We still want a solution of the form + +$$ +T(x,t) = f(x)g(t). +$$ + +We cannot use the boundary conditions at "infinity", because that would not help us in finding eigenvalues. Instead we can look at the functional description of the boundary condition at zero time. This is a Gaussian function in space and it is characterised by the width $L$. Hence the fact that the wire is infinitely long does not matter too much, because the temperature goes exponentially to zero for large distances. We can then choose a finite length $M$ and we will approximate the boundary condition with $2M$ as the length of the wire. + +The initial distribution is an even function in space. This means that an integral $-\infty0)$ or sink $(S<0)$, and substitution of Fourier's law in this equation gives + +$$ +\partial_t T - \kappa\nabla^2 T = S, +$$ (eq:Heat) + +where we have assumed constant heat conductivity. This is called a scalar diffusive field, or parabolic, equation, because in all source- and sink-free regions, a first-order time derivative of temperature is equal to a second-order spatial derivative of temperature. This is a general expression for a diffusion equation in three-dimensional space and time. A similar equation exists for other fields, for non-scalar fields as well, and we look at the electromagnetic field diffusion later. We first look at the heat equation in one and three dimensions. diff --git a/book/3_diffusion_fields/labs/figures/bar_boundary_conditions.svg b/book/3_diffusion_fields/labs/figures/bar_boundary_conditions.svg new file mode 100644 index 0000000..fed5034 --- /dev/null +++ b/book/3_diffusion_fields/labs/figures/bar_boundary_conditions.svg @@ -0,0 +1,1031 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/book/3_diffusion_fields/labs/figures/make_figures.py b/book/3_diffusion_fields/labs/figures/make_figures.py new file mode 100644 index 0000000..c4feb72 --- /dev/null +++ b/book/3_diffusion_fields/labs/figures/make_figures.py @@ -0,0 +1,43 @@ +"""Regenerate the schematic of the diffusive-fields lab. Run: python make_figures.py""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle, FancyArrowPatch + +INK, BAR, RES, OUT = "#6b7280", "#8a9199", "#00A6D6", "#EC6842" + + +def bar(ax, y, label): + ax.add_patch(Rectangle((0, y - 0.18), 1, 0.36, fc=BAR, ec=INK, alpha=0.35, lw=1)) + ax.text(0.5, y + 0.3, label, ha="center", va="bottom", color=INK, fontsize=10) + + +def arrow(ax, x0, x1, y): + ax.add_patch(FancyArrowPatch((x0, y), (x1, y), arrowstyle="-|>", mutation_scale=14, + color=OUT, lw=2)) + + +fig, ax = plt.subplots(figsize=(6.4, 3.0)) +# Dirichlet: ends clamped to reservoirs at the reference temperature +bar(ax, 1.25, "Example 1: $T(0,t)=T(L,t)=0$, ends held at the reference temperature") +for x0 in (-0.13, 1.0): + ax.add_patch(Rectangle((x0, 0.98), 0.13, 0.54, fc=RES, ec=RES, alpha=0.55)) +ax.text(-0.065, 0.9, "$T=0$", ha="center", va="top", color=INK, fontsize=9) +ax.text(1.065, 0.9, "$T=0$", ha="center", va="top", color=INK, fontsize=9) +arrow(ax, 0.1, -0.2, 1.25); arrow(ax, 0.9, 1.2, 1.25) +ax.text(1.24, 1.25, "heat can\nleave", ha="left", va="center", color=OUT, fontsize=8) + +# Neumann: insulated ends +bar(ax, 0.0, r"Example 2: $\partial_x T(0,t)=\partial_x T(L,t)=0$, insulated ends") +for x0 in (-0.05, 1.0): + ax.add_patch(Rectangle((x0, -0.27), 0.05, 0.54, fc="none", ec=INK, hatch="////", lw=1)) +ax.text(1.12, 0.0, "no flux", ha="left", va="center", color=INK, fontsize=8) + +ax.set_xticks([0, 1], ["$0$", "$L$"]); ax.set_yticks([]) +ax.set_xlim(-0.35, 1.45); ax.set_ylim(-0.5, 1.95) +for s in ("left", "right", "top"): + ax.spines[s].set_visible(False) +ax.spines["bottom"].set_color(INK); ax.tick_params(colors=INK) +ax.set_xlabel("$x$", color=INK) +fig.tight_layout() +fig.savefig("bar_boundary_conditions.svg", transparent=True, metadata={"Date": None}) diff --git a/book/3_diffusion_fields/labs/fwtools.py b/book/3_diffusion_fields/labs/fwtools.py new file mode 100644 index 0000000..e354ecf --- /dev/null +++ b/book/3_diffusion_fields/labs/fwtools.py @@ -0,0 +1,230 @@ +"""Plotting and self-check helpers for the ECTB2140 *Fields and Waves* labs, +diffusive-fields edition. + +You are **not** expected to read or edit this file during a lab session. It +exists so that your time goes on the physics -- series coefficients, decay +rates, Gaussian kernels -- rather than on plotting boilerplate. + +The self-check functions have the same names and behaviour as in the Lab 1 +and Lab 2 module. Everything here is one-dimensional in space: fields are +arrays of shape (nt, nx), one row per time. + +Requires: numpy, matplotlib, plotly. +""" + +from __future__ import annotations + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import TwoSlopeNorm +import plotly.graph_objects as go + +__all__ = [ + "show_profiles", "show_spacetime", "show_spacetime_slider", "show_records", + "check", "check_close", "check_abs", "check_scalar", +] + +# -------------------------------------------------------------------------- +# Figures +# -------------------------------------------------------------------------- + +def _fmt(v, fmt): + return fmt.format(v) if isinstance(fmt, str) else fmt(v) + + +def show_profiles(x, frames, values, *, value_name="t", value_fmt="{:.3g}", + unit="", reference=None, reference_name="initial condition", + xlabel="x [m]", ylabel="", title="", ylim=None, height=430, + ms_per_frame=90, colors=("#00A6D6", "#EC6842", "#7B61FF")): + """Profiles over time, with a play button and a draggable progress bar. + + ``frames`` is either one array of shape (len(values), nx), or a dict + ``{name: array}`` to draw several curves that evolve together. Press play + to run from the start, or drag the bar to any instant. Both run in the + browser; nothing here needs the kernel once the figure exists. + """ + if not isinstance(frames, dict): + frames = {"profile": frames} + frames = {k: np.asarray(v, float) for k, v in frames.items()} + x = np.asarray(x, float) + allv = np.concatenate([v.ravel() for v in frames.values()]) + lo, hi = (np.nanmin(allv), np.nanmax(allv)) if ylim is None else ylim + pad = 0.06 * (hi - lo if hi > lo else 1.0) + + fig = go.Figure() + if reference is not None: + fig.add_trace(go.Scatter(x=x, y=reference, mode="lines", name=reference_name, + line=dict(color="#8a9199", dash="dash", width=1.5))) + live = [] + for k, (name, arr) in enumerate(frames.items()): + fig.add_trace(go.Scatter(x=x, y=arr[0], mode="lines", name=name, + line=dict(color=colors[k % len(colors)], + width=2.5 if k == 0 else 2.0))) + live.append(len(fig.data) - 1) + + n_steps = len(values) + fig.frames = [go.Frame(name=str(k), + data=[go.Scatter(y=arr[k]) for arr in frames.values()], + traces=live) + for k in range(n_steps)] + still = dict(mode="immediate", frame=dict(duration=0, redraw=False), + transition=dict(duration=0)) + steps = [dict(method="animate", label=_fmt(v, value_fmt), args=[[str(k)], still]) + for k, v in enumerate(values)] + + fig.update_layout( + title=title, height=height, margin=dict(l=60, r=20, t=50, b=30), + xaxis_title=xlabel, yaxis_title=ylabel, + yaxis_range=[lo - pad, hi + pad], legend=dict(x=0.99, xanchor="right", y=0.99), + sliders=[dict(active=0, steps=steps, len=0.82, x=0.18, y=0, pad=dict(t=50, b=8), + currentvalue=dict(prefix=f"{value_name} = ", suffix=unit, + font=dict(size=13)))], + updatemenus=[dict(type="buttons", direction="left", showactive=False, + x=0.0, y=0, xanchor="left", yanchor="top", pad=dict(t=55, b=8), + buttons=[ + dict(label="▶ play", method="animate", + args=[None, dict(mode="immediate", fromcurrent=False, + frame=dict(duration=ms_per_frame, redraw=False), + transition=dict(duration=0))]), + dict(label="❚❚", method="animate", args=[[None], still])])]) + fig.show() + + +def show_spacetime(x, t, T, *, t_label="t [s]", x_label="x [m]", c_label="T [K]", + title="", signed=None, log_t=True, marks=None, curve=None, + curve_label="prediction", deeper_down=False): + """Colour map of T(x, t): time along the horizontal axis, position vertical. + + A signed field gets a diverging colour map centred on zero; a positive one + gets a sequential map. ``marks`` is an optional list of (t, x) points to + overlay, e.g. predicted peak times. + """ + T = np.asarray(T, float) + if signed is None: + signed = np.nanmin(T) < 0 < np.nanmax(T) + fig, ax = plt.subplots(figsize=(8.2, 3.6)) + if signed: + vmax = np.nanmax(np.abs(T)) + im = ax.pcolormesh(t, x, T.T, cmap="RdBu_r", shading="auto", + norm=TwoSlopeNorm(0.0, -vmax, vmax)) + else: + im = ax.pcolormesh(t, x, T.T, cmap="inferno", shading="auto", vmin=0) + if curve is not None: + ax.plot(curve[0], curve[1], "--", color="#00A6D6", lw=2, label=curve_label) + ax.legend(loc="upper left", fontsize=8) + if marks is not None: + mt, mx = zip(*marks) + ax.plot(mt, mx, "o", mfc="none", mec="#00A6D6", mew=1.5, ms=7, label="predicted peak") + ax.legend(loc="upper left", fontsize=8) + if log_t: + ax.set_xscale("log") + ax.set_xlim(np.min(t), np.max(t)) # an overlaid curve must not stretch the axis + if deeper_down: # depth axes read downward + ax.invert_yaxis() + ax.set_xlabel(t_label); ax.set_ylabel(x_label); ax.set_title(title) + fig.colorbar(im, ax=ax, label=c_label, pad=0.02) + fig.tight_layout() + plt.show() + + +def show_spacetime_slider(x, t, maps, values, *, value_name="a", value_fmt="{:.3g}", + unit="", t_label="t [s]", x_label="x [m]", c_label="T [K]", + title="", marks=None, height=430): + """One space-time map per slider position, for a parameter of the problem. + + ``maps`` has shape (len(values), len(t), len(x)). ``marks[k]`` optionally + draws a vertical line at that time on map k, e.g. a predicted merge time. + The slider updates the map and the line together, client-side. + """ + maps = np.asarray(maps, float) + x, t = np.asarray(x, float), np.asarray(t, float) + zmin, zmax = float(np.nanmin(maps)), float(np.nanmax(maps)) + + def line(k): + if marks is None or not np.isfinite(marks[k]) or marks[k] <= t[0]: + return [] + return [dict(type="line", x0=marks[k], x1=marks[k], y0=x[0], y1=x[-1], + xref="x", yref="y", line=dict(color="#00A6D6", width=2, dash="dash"))] + + fig = go.Figure(go.Heatmap(z=maps[0].T, x=t, y=x, colorscale="Inferno", + zmin=zmin, zmax=zmax, colorbar=dict(title=c_label))) + steps = [dict(method="update", label=_fmt(v, value_fmt), + args=[{"z": [m.T]}, {"shapes": line(k)}]) + for k, (v, m) in enumerate(zip(values, maps))] + fig.update_layout( + title=title, height=height, margin=dict(l=60, r=20, t=50, b=30), + xaxis=dict(title=t_label, type="log"), yaxis_title=x_label, + shapes=line(0), + sliders=[dict(active=0, steps=steps, len=0.82, x=0.12, y=0, pad=dict(t=50, b=8), + currentvalue=dict(prefix=f"{value_name} = ", suffix=unit, + font=dict(size=13)))]) + fig.show() + + +def show_records(t, records, *, t_label="t [s]", y_label="T [K]", title="", + peaks=None, log_t=True): + """Time series at fixed positions, as a thermometer or a receiver records them. + + ``records`` maps a label to an array over ``t``. ``peaks`` optionally maps + the same labels to a predicted (t_peak, value_peak), drawn as an open circle. + """ + fig, ax = plt.subplots(figsize=(8.2, 3.4)) + _ls = ("-", "--", "-.", ":") # not colour alone: the C2/C3 pair is a common confusion + for k, (lab, y) in enumerate(records.items()): + ax.plot(t, y, color=f"C{k}", ls=_ls[k % 4], lw=2, label=lab) + if peaks and lab in peaks: + ax.plot(*peaks[lab], "o", mfc="none", mec=f"C{k}", mew=1.5, ms=9) + if log_t: + ax.set_xscale("log") + ax.set_xlabel(t_label); ax.set_ylabel(y_label); ax.set_title(title) + ax.grid(alpha=0.3); ax.legend(fontsize=8) + fig.tight_layout() + plt.show() + + +# -------------------------------------------------------------------------- +# Self-checks (same behaviour as the Lab 1 and Lab 2 module) +# -------------------------------------------------------------------------- + +def check(label: str, ok: bool, hint: str = "") -> None: + """Report a pass, or raise with a hint about what to look at.""" + if ok: + print(f" [ok] {label}") + else: + raise AssertionError(f"{label} -- {hint}" if hint else label) + + +def check_close(label: str, got, want, rtol=0.05, where=None, hint: str = "") -> None: + """Compare two arrays where both are finite (and where `where` is True).""" + got, want = np.asarray(got, float), np.broadcast_to(np.asarray(want, float), np.shape(got)) + m = np.isfinite(got) & np.isfinite(want) + if where is not None: + m = m & where + if not m.any(): + raise AssertionError(f"{label} -- nothing left to compare; the mask removed every point") + rel = np.abs(got[m] - want[m]) / np.maximum(np.abs(want[m]), 1e-30) + worst = float(np.max(rel)) + check(f"{label}: worst error {worst:.2%}", worst < rtol, + hint or f"worst relative error {worst:.2%} exceeds {rtol:.0%}") + + +def check_abs(label: str, got, atol, where=None, hint: str = "") -> None: + """Compare an array against zero on an absolute scale, in the field's own units.""" + got = np.asarray(got, float) + m = np.isfinite(got) + if where is not None: + m = m & where + if not m.any(): + raise AssertionError(f"{label} -- nothing left to compare") + worst = float(np.max(np.abs(got[m]))) + check(f"{label}: worst |value| {worst:.2e}", worst < atol, + hint or f"worst deviation {worst:.2e} exceeds {atol:.1e}") + + +def check_scalar(label: str, got: float, want: float, rtol: float = 0.01, + unit: str = "", hint: str = "") -> None: + """Compare two single numbers and report the relative discrepancy.""" + got, want = float(got), float(want) + rel = abs(got - want) / max(abs(want), 1e-30) + check(f"{label}: {got:.4g}{unit} vs {want:.4g}{unit} ({rel:.2%} apart)", + rel < rtol, hint or f"these should agree to better than {rtol:.0%}") diff --git a/book/3_diffusion_fields/labs/week05_diffusion.md b/book/3_diffusion_fields/labs/week05_diffusion.md new file mode 100644 index 0000000..52a6600 --- /dev/null +++ b/book/3_diffusion_fields/labs/week05_diffusion.md @@ -0,0 +1,788 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +mystnb: + # Workbook page: the task cells contain `___` blanks by design, so it must + # not be executed at build time. Readers run it themselves with Live Code. + execution_mode: 'off' +--- + +# Lab 5: Diffusive Fields + +:::{admonition} Computer lab +:class: note + +A practical companion to the lecture on [the heat equation in one dimension](../heat_equation_1d.md). Each task states a physical question, gives the steps, and ends with a self-check you can run. + +Everything you write here is plain NumPy. The module `fwtools` only draws and checks: `fw.show_profiles` (a profile per instant, with a play button), `fw.show_spacetime` (a colour map over position and time), `fw.show_records` (time series at fixed positions) and the `fw.check...` functions. Every slider runs in the browser. +::: + +## Learning objectives + +- **Decaying modes.** Build an initial temperature from a Fourier series and explain why mode $n$ decays at $\kappa(n\pi/L)^2$, so fine structure goes first. +- **Boundaries.** Ends held at a fixed temperature let heat leave; insulated ends conserve it and drive the bar to its mean. +- **Superposition.** Evolve any initial temperature with the Green's function, and turn the $t\propto x^2$ it implies into a measurement: of $\kappa$ in a bar, of depth in the Earth. + +--- + +## Part 0 — Setup + +Run this once. It contains no physics: it fetches two packages the browser lacks and locates `fwtools`. + +```{code-cell} ipython3 +import sys, pathlib + +import numpy as np +import matplotlib.pyplot as plt +from scipy.constants import mu_0, epsilon_0 + +# --- Live Code housekeeping, not part of the physics ----------------------- +try: + import plotly.io as pio +except ModuleNotFoundError: + print("Fetching plotly. A few seconds, and only the first time...") + import micropip + await micropip.install("plotly") + import plotly.io as pio + +try: + import nbformat # noqa: F401 +except ModuleNotFoundError: + import micropip, types + try: + await micropip.install("nbformat") + except Exception: + _nb = types.ModuleType("nbformat") + _nb.__version__ = "5.10.4" + sys.modules["nbformat"] = _nb + +for _p in (".", "book/3_diffusion_fields/labs"): + if (pathlib.Path(_p) / "fwtools.py").exists(): + sys.path.insert(0, _p) + break +try: + import fwtools as fw +except ModuleNotFoundError: + from pyodide.http import pyfetch + _r = await pyfetch("fwtools.py") + pathlib.Path("fwtools.py").write_bytes(await _r.bytes()) + import fwtools as fw + +pio.renderers.default = "plotly_mimetype+notebook" +# --------------------------------------------------------------------------- +print("ready") +``` + +--- + +## Part 1 — A bar of finite length + +A copper bar of length $L = 0.5$ m has thermal diffusivity $\kappa = 1.1\times10^{-4}$ m²/s. Its temperature $T$ is measured from a reference value, so it can be negative. Without sources, the temperature obeys + +$$ \partial_t T = \kappa\,\partial_x^2 T . $$ + +The only time scale the bar has is $L^2/\kappa \approx 38$ min, so the natural clock is the Fourier number $\mathrm{Fo} = \kappa t/L^2$. The two ways of closing the ends lead to different series and to a different final state. + +```{figure} figures/bar_boundary_conditions.svg +:name: fig-lab5-bcs +:width: 80% + +The two bars of this part. In Tasks 1 and 2 the ends are clamped to reservoirs at the reference temperature. In Task 3 they are insulated. +``` + +### Task 1 — build the initial temperature from sines + +With both ends held at $T=0$, separation of variables gives + +$$ T(x,t) = \sum_{n=1}^{\infty} C_n \sin\!\left(\frac{n\pi x}{L}\right) \exp\!\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], \qquad C_n = \frac{2}{L}\int_0^L h(x)\sin\!\left(\frac{n\pi x}{L}\right)\mathrm{d}x , $$ + +where the coefficient formula follows from the orthogonality of the sines. At $t=0$ the series must rebuild the initial temperature $h(x)$. Here $h$ is $0$, $-1$, $1$, $0$ K on the four quarters of the bar. + +1. Compute all $C_n$ at once by projecting $h$ onto every row of `S`. +2. Write the partial sum of the first `m` modes. +3. Drag the slider. **Before you do**, predict which coefficients vanish, using the symmetry of $h$ about $x=L/2$. + +```{code-cell} ipython3 +L, kappa = 0.5, 1.1e-4 # copper bar: length [m], thermal diffusivity [m^2/s] +ncell, N = 1200, 400 # grid cells, number of modes kept +dx = L / ncell +x = (np.arange(ncell) + 0.5) * dx # cell centres, so an integral is a plain sum times dx +n = np.arange(1, N + 1) +S = np.sin(np.outer(n, np.pi * x / L)) # row n-1 is sin(n pi x / L), shape (N, ncell) +h = np.where(x < L/4, 0.0, np.where(x < L/2, -1.0, + np.where(x < 3*L/4, 1.0, 0.0))) # initial temperature [K] +print(f"L^2/kappa = {L**2/kappa:.0f} s; slowest possible mode, tau_1 = L^2/(kappa pi^2)" + f" = {L**2/(kappa*np.pi**2):.0f} s") + +# Task 1 -- two blanks. An integral over the bar is np.sum(...) * dx; summing along +# axis=1 of an (N, ncell) array does all N integrals at once. +C = ___ # all C_n, shape (N,) +partial = lambda m: ___ # sum of the first m modes at t = 0, shape (ncell,) + +# --- given: the partial sums, one per slider position --- +Ns = np.array([1, 2, 4, 6, 10, 20, 40, 100, 400]) +fw.show_profiles(x, np.array([partial(m) for m in Ns]), Ns, value_name="modes kept", + value_fmt="{:d}", reference=h, reference_name="h(x)", + ylabel="T [K]", title="Partial sums of the sine series") +print("C_1 .. C_6 =", np.round(C[:6], 4)) + +# --- self-check (leave this alone) --- +_C_exact = 2 * (2*np.cos(n*np.pi/2) - np.cos(n*np.pi/4) - np.cos(3*n*np.pi/4)) / (np.pi * n) +fw.check_abs("the first 40 C_n match the closed form", C[:40] - _C_exact[:40], 1e-3, + hint="C_n = (2/L) * integral of h(x) sin(n pi x/L) dx, for every n at once") +_away = np.min(np.abs(x[:, None] - np.array([L/4, L/2, 3*L/4])), axis=1) > 0.04 +fw.check_abs("400 modes rebuild h(x) more than 4 cm from a jump", partial(N) - h, 0.02, where=_away, + hint="partial(m) should weight the first m rows of S by the first m coefficients") +``` + +:::{admonition} Solution — Task 1 +:class: dropdown + +```python +C = (2 / L) * np.sum(h * S, axis=1) * dx +partial = lambda m: C[:m] @ S[:m] +``` +::: + +:::{admonition} Only even modes appear, and the overshoot at a jump never shrinks +:class: important dropdown + +$h$ is antisymmetric about $x=L/2$, while $\sin(n\pi x/L)$ is symmetric about that point for odd $n$ and antisymmetric for even $n$. The projection of an antisymmetric function onto a symmetric one is zero, so $C_1=C_3=C_5=\dots=0$. + +Near each jump the partial sums overshoot by about 9% of the jump. At $x=L/2$ the jump is 2 K, so the peak reaches about 1.18 K, whether you keep 10 modes or 100. Adding terms squeezes the overshoot towards the jump without lowering it; at 400 modes it becomes narrower than the grid spacing, and the plotted peak reads slightly lower. This is the Gibbs phenomenon: the series converges at every point, but not uniformly near a discontinuity. +::: + +### Task 2 — let it diffuse + +Every mode keeps its shape and decays with its own time constant + +$$ \tau_n = \frac{L^2}{\kappa\,n^2\pi^2} . $$ + +Mode $n$ dies $n^2$ times faster than mode 1. **Predict** which mode is the last one visible, given the coefficients of Task 1. + +1. Build the decay factors for all times and all modes in one array. +2. Weight each mode by $C_n$ and by its decay, and sum the modes. + +```{code-cell} ipython3 +t = np.logspace(-1, np.log10(0.25 * L**2 / kappa), 90) # 0.1 s .. 568 s, i.e. Fo = 0.25 + +# Task 2 -- two blanks. +D = ___ # exp(-kappa (n pi/L)^2 t), shape (len(t), N) +T = ___ # T(x, t), shape (len(t), nx) + +# --- given: the space-time map, then the profiles as a movie --- +fw.show_spacetime(x, t, T, t_label="t [s]", x_label="x [m]", c_label="T [K]", + title="Ends held at T = 0") +fw.show_profiles(x, T[::2], t[::2], value_name="t", value_fmt="{:.3g}", unit=" s", + reference=h, ylabel="T [K]", title="Ends held at T = 0", ylim=(-1.2, 1.2)) + +# --- self-check (leave this alone) --- +fw.check("T has one row per time and one column per grid cell", np.shape(T) == (len(t), ncell)) +_i1, _i2 = np.argmin(np.abs(t - 0.05*L**2/kappa)), np.argmin(np.abs(t - 0.10*L**2/kappa)) +_A = np.max(np.abs(T), axis=1) +_rate = np.log(_A[_i1] / _A[_i2]) / (t[_i2] - t[_i1]) +fw.check_scalar("late-time decay rate of your T equals kappa (2 pi/L)^2", _rate, + kappa * (2*np.pi/L)**2, rtol=0.01, unit=" 1/s", + hint="each mode needs its own exponent kappa (n pi/L)^2 t") +``` + +:::{admonition} Solution — Task 2 +:class: dropdown + +```python +D = np.exp(-kappa * (n * np.pi / L)**2 * t[:, None]) +T = (C * D) @ S +``` +::: + +:::{admonition} Mode 2 survives, because mode 1 was never there +:class: important dropdown + +The sharp corners of $h$ are made of high modes, and those vanish within seconds: after one second, modes above $n\approx 15$ have decayed by a factor of $e$. By $t\approx 60$ s only $n=2$ and a trace of $n=4$ remain, so the profile has become a single sine wave that shrinks without changing shape. Mode 1 would have decayed four times more slowly still, but $h$ contains none of it. The heat flows out through the clamped ends and between the warm and cold halves, and the bar returns to the reference temperature. + +The measured decay rate is the second eigenvalue, $\kappa(2\pi/L)^2 = 0.0174$ s$^{-1}$. A long-time temperature record thus reveals the lowest mode present, and with it $\kappa/L^2$. +::: + +### Task 3 — insulate the ends + +With $\partial_xT = 0$ at both ends, no heat crosses them. The eigenfunctions become cosines, and $\lambda=0$ is now an eigenvalue: + +$$ T(x,t) = C_0 + \sum_{n=1}^{\infty} C_n \cos\!\left(\frac{n\pi x}{L}\right)\exp\!\left[-\kappa\left(\frac{n\pi}{L}\right)^2 t\right], \qquad C_0 = \frac{1}{L}\int_0^L h\,\mathrm{d}x, \quad C_n = \frac{2}{L}\int_0^L h\cos\!\left(\frac{n\pi x}{L}\right)\mathrm{d}x . $$ + +The initial temperature is now $1$, $-1$, $1$ K on the three thirds of the bar. The cosine projection is given; write $C_0$, the one mode the clamped bar of Task 2 did not have. **Predict** the temperature at late times, and whether the heat content $\int_0^L T\,\mathrm{d}x$ changes. + +```{code-cell} ipython3 +Cc = np.cos(np.outer(n, np.pi * x / L)) # row n-1 is cos(n pi x / L) +h_ins = np.where(x < L/3, 1.0, np.where(x < 2*L/3, -1.0, 1.0)) # initial temperature [K] + +# Task 3 -- one blank: the n = 0 coefficient, the mode the clamped bar did not have. +C0 = ___ # the n = 0 coefficient [K] + +# --- given: the rest is Task 2 with cosines in place of sines --- +Cn = (2 / L) * np.sum(h_ins * Cc, axis=1) * dx # C_n for n = 1..N, shape (N,) +T_ins = C0 + (Cn * D) @ Cc # T(x, t), reusing D from Task 2 + +Q = np.sum(T_ins, axis=1) * dx # heat content, per unit area and heat capacity [K m] + +# --- given: the same two views as Task 2, to be compared with it --- +fw.show_spacetime(x, t, T_ins, t_label="t [s]", x_label="x [m]", c_label="T [K]", + title="Insulated ends") +fw.show_profiles(x, T_ins[::2], t[::2], value_name="t", value_fmt="{:.3g}", unit=" s", + reference=h_ins, ylabel="T [K]", title="Insulated ends", ylim=(-1.2, 1.2)) +print(f"heat content: min {Q.min():.6f}, max {Q.max():.6f} K m over {t[0]:.1f} s .. {t[-1]:.0f} s") + +# --- self-check (leave this alone) --- +fw.check_close("heat content stays at C0 L = 1/6 K m at every time", Q, (1/3) * L, rtol=1e-3) +_away3 = np.min(np.abs(x[:, None] - np.array([L/3, 2*L/3])), axis=1) > 0.02 +fw.check_abs("at t = 0.1 s the series still equals h(x) away from the jumps", + T_ins[0] - h_ins, 0.01, where=_away3, + hint="compare your cosine projection with the sine projection of Task 1") +``` + +:::{admonition} Solution — Task 3 +:class: dropdown + +```python +C0 = np.sum(h_ins) * dx / L +``` +::: + +:::{admonition} Insulated ends conserve the heat, and the bar settles at its mean +:class: important dropdown + +Every cosine with $n\ge1$ integrates to zero over the bar, so the heat content is $C_0L$ at all times, fixed by $h$ alone. The mode $n=0$ has decay rate zero, so it is all that remains: $T\to C_0 = 1/3$ K everywhere. In Task 2 every eigenvalue was positive and the temperature went to zero, because the clamped ends could absorb or supply heat. The boundary condition, not the equation, decides where the heat ends up. +::: + +--- + +## Part 2 — An infinitely long wire + +On an infinite wire the lecture notes solve the same equation twice. First in closed form, for an initial Gaussian $T(x,0) = T_0\exp(-x^2/L^2)$: + +$$ T(x,t) = \frac{T_0 L}{\sqrt{4\kappa t + L^2}}\exp\!\left(-\frac{x^2}{4\kappa t + L^2}\right). $$ + +Then again as a cosine series, by putting the wire in a box of length $2M$ and letting $M$ grow. Task 7 takes up that second route; Parts 2 and 3 use the closed form, because it is the one you can read. + +In this part and again in Task 7, $L$ is the width of the initial Gaussian and not a length of wire; the lecture notes use the one letter for both. The code keeps them apart: `L` is the bar of Part 1, `Lg` the Gaussian here, `Lw` the Gaussian of Task 7. Time enters only through $\kappa t$, which has units of m², so every result below holds for any material. The figures below use $L=3$ m, $-10 kt[0]}, + t_label=r"$\kappa t$ [m$^2$]", y_label="T [K]", + title="Thermometer records; circles mark the predicted peaks") + +# --- self-check (leave this alone) --- +_xw = np.linspace(-2000.0, 2000.0, 200001) +_heat = [np.sum(T_gauss(_xw, k)) * (_xw[1] - _xw[0]) for k in (0.1, 10.0, 1000.0)] +fw.check_close("heat content T0 L sqrt(pi) at kappa t = 0.1, 10, 1000 m^2", + _heat, T0 * Lg * np.sqrt(np.pi), rtol=1e-3, + hint="the prefactor must shrink exactly as fast as the width grows") +_ktf = np.logspace(-1, 3, 4001) +_rec = T_gauss(5.0, _ktf) +fw.check_scalar("your kt_peak(5 m) against the peak of your own record", kt_peak(5.0), + _ktf[np.argmax(_rec)], rtol=0.01, unit=" m^2") +fw.check_scalar("peak temperature at 5 m", np.max(_rec), T0 * Lg / (np.sqrt(2*np.e) * 5.0), + rtol=0.01, unit=" K") +``` + +:::{admonition} Solution — Task 4 +:class: dropdown + +```python +def T_gauss(x, kt): + return T0 * Lg / np.sqrt(4*kt + Lg**2) * np.exp(-x**2 / (4*kt + Lg**2)) + +def kt_peak(x): + return (2*x**2 - Lg**2) / 4 +``` +::: + +:::{admonition} The heat spreads without loss, and a farther thermometer answers later and weaker +:class: important dropdown + +The width grows as $\sqrt{4\kappa t+L^2}$, and the amplitude falls by exactly the same factor, so the area under the profile stays $T_0L\sqrt\pi$. With no ends, no heat is lost; it is only redistributed. At late times the temperature everywhere falls as $t^{-1/2}$. + +The thermometers show the two faces of diffusion. The peak arrives at $\kappa t_p\propto x^2$: twice as far away takes four times as long. Its height falls as $1/\lvert x\rvert$. Inside $\lvert x\rvert2a^2$ the profile has a single maximum, and it becomes ever harder to distinguish from one Gaussian of twice the heat. Sources that are closer together merge sooner, as $a^2$. For $a thresh)]) +t_dep = np.array(t_dep) + +slope, intercept = np.polyfit(xcs**2, t_dep, 1) +kappa_fit = 1 / (4 * slope * np.log(1 / thresh)) +plt.figure(figsize=(5.2, 3.2)) +plt.plot(xcs**2, t_dep, "o", color="#00A6D6", label="measured") +plt.plot(xcs**2, slope * xcs**2 + intercept, "-", color="#8a9199", + label=f"fit: $\\kappa$ = {kappa_fit:.2e} m$^2$/s") +plt.xlabel("$x_c^2$ [m$^2$]"); plt.ylabel("departure time [s]") +plt.grid(alpha=0.3); plt.legend(fontsize=8); plt.tight_layout(); plt.show() +print(f"departure times [s]: {np.round(t_dep, 2)}") +print(f"intercept: {intercept:.2f} s (-w^2/(4 kappa) = {-w_b**2/(4*kappa):.2f} s)") +print(f"kappa from the fit: {kappa_fit:.3e} m^2/s (the bar was built with {kappa:.3e})") + +# --- self-check (leave this alone) --- +_h_check = np.exp(-(x - 0.15)**2 / w_b**2)[::4] +fw.check_abs("at t = 0.1 s the superposition reproduces the series solution", + T_free(t6[0], _h_check) - bar[0], 1e-6, + hint="the kernel must carry unit heat: check the 1/sqrt(4 pi kappa t) factor " + "and that the sum over x' is multiplied by dxs") +fw.check_scalar("kappa recovered from the departure times", kappa_fit, kappa, rtol=0.05, + unit=" m^2/s") +``` + +:::{admonition} Solution — Task 6 +:class: dropdown + +```python +def T_free(tv, h_on_grid): + G = np.exp(-d2 / (4 * kappa * tv)) / np.sqrt(4 * np.pi * kappa * tv) + return G @ h_on_grid * dxs +``` +::: + +:::{admonition} The ends are felt when the tail reaches them, and that measures $\kappa$ +:class: important dropdown + +The bump spreads to width $\sqrt{4\kappa t + w^2}$, so the temperature at the near end, as a +fraction of the peak, is $\exp[-x_c^2/(4\kappa t + w^2)]$. Setting that ratio to the threshold +$\theta$ gives + +$$ 4\kappa\,t_{\text{dep}} = \frac{x_c^2}{\ln(1/\theta)} - w^2 , $$ + +a straight line in $x_c^2$, which is what the second figure shows. Its slope is +$1/(4\kappa\ln(1/\theta))$, so the fit returns $\kappa$ to about 1%. The intercept is negative: +$-w^2/(4\kappa) = -2.05$ s predicted against $-1.96$ s measured, the head start the bump gets from +the width it already had at $t=0$. + +$\theta$ cancels out of $\kappa$ because the Gaussian tail is exactly log-quadratic, so moving the +threshold scales every departure time by one common factor and the slope absorbs it. Set `thresh` +to 0.002 and to 0.02: the nearest bump departs at 1.7 s instead of 4.0 s, yet $\kappa$ comes back as +$1.08\times10^{-4}$ and $1.09\times10^{-4}$ m²/s against the $1.09\times10^{-4}$ of the 1% run. A bump +whose tail was not Gaussian would not cancel. Push `thresh` to 0.05 and the fit finally drifts 5% +low, because by then the departure is being read off the shoulder of the bump rather than its tail. + +Two things are worth taking from this. A boundary makes itself felt through the tail of the +distribution, long before the diffusion length $\sqrt{4\kappa t}$ equals the distance to it. +And the same $t\propto x^2$ that set the peak time in Task 4 sets this departure time, because +a diffusive field has only one way to convert a distance into a time. +::: + + +### Task 7 — how big a box, and how many modes? + +Part 2 used the closed form on the infinite wire. The lecture notes reach the same field a second +way: put the wire in a box of length $2M$, expand the initial Gaussian in cosines, and let each +mode decay at its own rate. With $M\gg L$ the coefficients can be integrated in closed form, and +the series collapses to + +$$ T(x,t) = \frac{T_0L\sqrt{\pi}}{M}\left\{\frac{1}{2} + \sum_{n=1}^{N}\cos\!\left(\frac{n\pi x}{M}\right)\exp\!\left[-\left(\frac{n\pi}{2M}\right)^2\!(4\kappa t+L^2)\right]\right\} . $$ + +The box is fictitious, so $M$ and $N$ are ours to choose, and a wrong choice shows up as an error +against the exact solution. Take $L=10$ m, $-3L11}{'tau_r [s]':>11}{'t_p / tau_r':>13}") +for k, sv in media.items(): + print(f"{k:<30}{t_peak_em(z_s, sv):>11.3g}{eps_r*epsilon_0/sv:>11.3g}" + f"{peak_over_relax(z_s, sv):>13.3g}") + +print("\na survey looking 500 m down, able to record for 10 ms:") +for k, sv in media.items(): + _tp = t_peak_em(500.0, sv) + print(f" {k:<30} t_p = {_tp:.3g} s {'inside' if _tp < 1e-2 else 'beyond'} the window") + +# --- self-check (leave this alone) --- +_tf = np.logspace(-8, 1, 20001) +for k, sv in media.items(): + fw.check_scalar(f"t_p at {z_s:.0f} m in {k.split(',')[0]}", t_peak_em(z_s, sv), + _tf[np.argmax(np.abs(E_step(z_s, _tf, sv)))], rtol=0.01, unit=" s") +fw.check("basement to seawater is a factor 4e4 in arrival time", + np.isclose(t_peak_em(z_s, 4.0) / t_peak_em(z_s, 1e-4), 4e4, rtol=1e-9), + "t_p is linear in sigma, and the two conductivities differ by 4e4") +``` + +:::{admonition} Solution — Task 8 +:class: dropdown + +```python +def t_peak_em(z, s=sigma): + return s * mu_0 * z**2 / 2 + +def peak_over_relax(z, s): + return s**2 * mu_0 * z**2 / (2 * eps_r * epsilon_0) +``` +::: + +:::{admonition} Depth maps onto arrival time, as $t_p\propto z^2$ +:class: important dropdown + +A layer ten times deeper answers a hundred times later: 63 µs at 100 m, 6.3 ms at 1 km, for $\sigma=0.01$ S/m. On the map your curve $t_p(z)$ runs along the shoulder of the diffusing field, because $|E_x|$ has fallen to $1/\sqrt e$ of its surface value exactly where $t = \sigma\mu z^2/2$. Time-domain electromagnetic soundings rest on this: the later part of a record senses the deeper Earth. In a more conductive Earth the arrival is later, since $t_p\propto\sigma$. The ratio $t_p/\tau_r$ is 7000 at 100 m and grows as $z^2$, so neglecting $\varepsilon\,\partial_t\boldsymbol{E}$ is well justified at these depths. The blurring of Task 5 applies here too: deep structure arrives late, and therefore blurred. + +The peak amplitude is 0.24 V/m at 100 m and 0.024 V/m at 1 km. A factor of ten in depth costs a factor of ten in signal and a factor of a hundred in time, which is the $1/\lvert x\rvert$ law of Task 4 with depth in place of distance along the wire. + +The chapter's three-dimensional impulse response peaks at $t_p=\sigma\mu r^2/6$, three times earlier than the $\sigma\mu z^2/2$ you found here, and its amplitude falls as $1/r^3$. Both follow from the same differentiation: a prefactor $t^{-p}$ gives $t_p=\sigma\mu r^2/(4p)$, with $p=1/2$ for the sheet and $p=3/2$ for the point. +::: + +:::{admonition} Conductivity sets the clock, and it sets the margin twice over +:class: important dropdown + +At 300 m the peak arrives after 0.23 s in seawater, 0.57 ms in sediments and 5.7 µs in basement. That is a factor $4\times10^{4}$ from basement to seawater, exactly the ratio of the two conductivities, because $t_p=\sigma\mu z^2/2$ is linear in $\sigma$. A conductive ground holds the field up; a resistive one lets it through almost at once. + +The margin behaves differently. $\tau_r=\varepsilon/\sigma$ grows as the ground becomes resistive while $t_p$ shrinks, so + +$$ \frac{t_p}{\tau_r} = \frac{\sigma^2\mu z^2}{2\varepsilon} $$ + +carries $\sigma$ twice and falls a hundredfold for every tenfold drop in conductivity: $10^{10}$ in seawater, $6.4\times10^{4}$ in sediments, and **6.4 in basement**. The diffusive approximation is comfortable in the first two and running out in the third, and it is the *resistive* ground that breaks it, not the conductive one. Dropping $\varepsilon\,\partial_t\boldsymbol{E}$ was never a statement about depth. It is a statement about how fast charge can relax, and in a resistive ground it cannot relax quickly. + +The survey question resolves the same way from the other end. Looking 500 m down with a 10 ms record: seawater peaks at 0.63 s, far beyond the window, so the target is never reached in time. Sediments peak at 1.6 ms, comfortably inside. Basement peaks at 16 µs, inside the window but so early that a real transmitter is still switching off, and the useful signal is buried in the turn-off transient. Conductivity, not instrument quality, decides which depths a time-domain sounding can reach. +::: diff --git a/book/_toc.yml b/book/_toc.yml index c509667..43b5ab1 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -96,6 +96,13 @@ parts: # - file: extension_exercises/014.md # - file: extension_exercises/015.md # - file: extension_exercises/016.md + - caption: Diffusive fields + chapters: + - file: 3_diffusion_fields/intro.md + sections: + - file: 3_diffusion_fields/heat_equation_1d.md + - file: 3_diffusion_fields/em_diffusion.md + - file: 3_diffusion_fields/labs/week05_diffusion.md - caption: Miscellaneous chapters: - file: references.md