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

178 5 SolvingPartialDifferentialEquations Exercise5.2:Compute temperaturevariations in theground Thesurface temperatureat thegroundshowsdailyandseasonaloscillations.When the temperature rises at the surface, heat is propagated into the ground, and the coefficientˇ in the diffusion equation determines how fast this propagation is. It takes some time before the temperature rises down in the ground. At the surface, the temperature has then fallen. We are interested in how the temperature varies downin thegroundbecauseof temperatureoscillationson thesurface. Assuminghomogeneoushorizontalpropertiesof theground,at least locally,and novariations of the temperature at the surface at a fixedpoint of time,we canne- glect thehorizontalvariationsof the temperature.Thenaone-dimensionaldiffusion equation governs the heat propagation along a vertical axis called x. The surface corresponds tox D 0 and thex axis point downwards into the ground. There is nosource termin theequation(actually, if rocks in thegroundare radioactive, they emit heat and that can bemodeled by a source term, but this effect is neglected here). At somedepthxDLweassumethat theheatchanges inxvanish, so@u=@xD 0 is an appropriate boundarycondition atx D L. Weassume a simple sinusoidal temperaturevariationat thesurface: u.0;t/DT0CTasin 2 P t ; whereP is theperiod, takenhereas24hours (24 60 60s). Theˇ coefficientmay be set to10 6m2=s. Time is thenmeasured in seconds. Set appropriatevalues for T0 andTa. a) Showthat thepresentproblemhasananalytical solutionof the form u.x;t/DACBe rx sin.!t rx/; forappropriatevaluesofA,B,r, and!. b) Solve thisheatpropagationproblemnumerically forsomedaysandanimate the temperature. Youmay use the Forward Euler method in time. Plot both the numerical and analytical solution. As initial condition for the numerical solu- tion, use the exact solution duringprogramdevelopment, andwhen the curves coincide in theanimationforall times,your implementationworks,andyoucan then switch to a constant initial condition: u.x;0/ D T0. For this latter initial condition,howmanyperiodsofoscillationsarenecessarybefore there isagood (visual)matchbetween thenumerical andexact solution (despitedifferencesat t D0)? Filename:ground_temp.py. Exercise5.3:Compare implicitmethods Anequally stable, butmore accuratemethod than theBackwardEuler scheme, is the so-called 2-step backward scheme, which for an ODE u0 D f.u;t/ can be expressedby 3unC1 4unCun 1 2 t Df.unC1;tnC1/: