Page - 307 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

9.3 Exercises 307 Theinitialconditionis thefamousandwidelyusedGaussianfunctionwithstandard deviation (or “width”)σ, which is here taken to be small,σ = 0.01, such that the initial condition is a peak.Thispeakwill thendiffuseandbecomelowerandwider. Computeu(x,t)untilubecomesapproximatelyconstantover the domain. Filename:gaussian_diffusion.py. Remarks Running the simulation with σ = 0.2 results in a constant solution u ≈ 1 as t → ∞, while one might expect from “physics of diffusion” that the solution should approach zero. The reason is that we apply Neumann conditions as boundary conditions. One can then easily show that the area under theu curve remainsconstant. Integrating thePDE gives ∫ 1 −1 ∂u ∂t dx=β ∫ 1 −1 ∂2u ∂x2 dx. Using the Gauss divergence theorem on the integral on the right-hand and moving the time-derivativeoutside the integralon the left-handside results in ∂ ∂t ∫ 1 −1 u(x,t)dx=β [ ∂u ∂x ]1 −1 =0. (Recall that∂u/∂x = 0 at the end points.) The result means that∫1−1udx remains constant during the simulation. Giving the PDE an interpretation in terms of heat conduction can easily explain the result: with Neumann conditions no heat can escapefromthedomainsotheinitialheatwill justbeevenlydistributed,butnotleak out, so the temperature cannot go to zero (or the scaled and translated temperature u, to beprecise).Theareaunder the initial condition is 1, sowith asufficientlyfine mesh,u→1, regardlessofσ. Exercise9.7:VectorizeaFunctionforComputing theAreaofaPolygon Vectorize the implementation of the function for computing the area of a polygon in Exercise 5.6. Make a test function that compares the scalar implementation in Exercise 5.6 and the new vectorized implementation for the test cases used in Exercise5.6. Hint Notice that the formula x1y2 + x2y3 + ··· + xn−1yn=∑n−1i=0 xiyi+1 is the dot product of two vectors, x[:-1] and y[1:], which can be computed as numpy.dot(x[:-1], y[1:]), or more explicitly as numpy.sum(x[:-1] *y[1:]). Filename:polyarea_vec.py. Exercise9.8:ExploreSymmetry One can observe (and also mathematically prove) that the solution u(x,t) of the problem in Exercise 9.6 is symmetric aroundx = 0:u(−x,t) = u(x,t). In such a case, we can split the domain in two and compute u in only one half, [−1,0] or [0,1]. At the symmetry line x = 0 we have the symmetry boundary condition