Page - 181 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

5.2 Exercises 181 Remarks Runningthesimulationwith D0:2 results inaconstantsolutionu 1 as t ! 1, while onemight expect from “physics of diffusion” that the solution should approachzero. The reason is thatweapplyNeumannconditionsas bound- ary conditions. One can then easily show that the area under theu curve remains constant. Integrating thePDEgives 1Z 1 @u @t dxDˇ 1Z 1 @d2u @x2 dx: Using theGauss divergence theoremon the integral on the right-handandmoving the time-derivativeoutside the integralon the left-handside results in @ @t 1Z 1 u.x;t/dxDˇ @du @x 1 1 D0: (Recall that@u=@x D 0 at the endpoints.) The resultmeans thatR1 1udx remains constant during the simulation. Giving the PDEan interpretation in terms of heat conduction can easily explain the result: with Neumann conditions no heat can escapefromthedomainsotheinitialheatwill justbeevenlydistributed,butnotleak out, so the temperature cannotgo to zero (or the scaled and translated temperature u, tobeprecise). Theareaunder the initial condition is1, sowithasufficientlyfine mesh,u!1, regardlessof . Exercise5.7:Vectorizea function forcomputing theareaofapolygon Vectorize the implementation of the function for computing the area of a polygon in Exercise 2.5. Make a test function that compares the scalar implementation in Exercise2.5and thenewvectorized implementationfor the test casesused inExer- cise2.5. Hint Notice that the formula x1y2 C x2y3 C C xn 1yn D Pn 1iD0xiyiC1 is the dot product of two vectors, x[:-1] and y[1:], which can be computed as numpy.dot(x[:-1], y[1:]),ormoreexplicitlyasnumpy.sum(x[:-1]*y[1:]). Filename:polyarea_vec.py. Exercise5.8:Exploresymmetry One can observe (and alsomathematically prove) that the solutionu.x;t/ of the problem inExercise 5.6 is symmetric aroundx D 0: u. x;t/ Du.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 linex D 0wehave the symmetry boundary condition @u=@x D 0. Reformulate the problem inExercise 5.6 such thatwe compute only forx 2 Œ0;1 . Display the solution andobserve that it equals the right part of the solution inExercise5.6. Filename:symmetric_gaussian_diffusion.py.