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

172 5 SolvingPartialDifferentialEquations Scalingmeans that we introduce dimensionless independent and dependent variables, heredenotedbyabar: NuD u u uc u ; NxD x xc ; Nt D t tc ; whereuc isacharacteristicsizeof the temperature,u is somereferencetemper- ature,whilexc and tc arecharacteristic timeandspace scales. Here, it isnatural to chooseu as the initial condition, and setuc to the stationary (end) temper- ature. Then Nu 2 Œ0;1 , starting at 0 and ending at 1 as t ! 1. The lengthL isxc,while choosing tc ismore challenging, but one can argue for tc DL2=ˇ. The resultingequation for Nu reads @Nu @Nt D @2 Nu @Nx2; Nx2 .0;1/: Note that in this equation, thereareno physical parameters! Inotherwords,we have foundamodel that is independentof the lengthof the rodand thematerial it ismadeof (!). Wecaneasily solve this equationwithourprogrambysettingˇD1,LD1, I.x/ D 0, and s.t/ D 1. It turns out that the total simulation time (to “infin- ity”) can be taken as 1.2. Whenwehave the solution Nu.Nx; Nt/, the solutionwith dimensionKelvin, reflecting the true temperature inourmedium, isgivenby u.x;t/Du C.uc u /Nu.x=L;tˇ=L2/: Through this formulawe can quickly generate the solutions for a rodmade of aluminum,wood,or rubber - it is just amatterofplugging in the rightˇvalue. Figure5.3showsfoursnapshotsof thescaled(dimensionless) solution N.Nx; Nt/. Thepowerofscalingis toreducethenumberofphysicalparametersinaprob- lem,and in thepresentcase,we foundonesingleproblemthat is independentof thematerial (ˇ) and thegeometry (L). 5.1.5 Vectorization Occasionally in this book,we showhow to speedupcodeby replacing loopsover arraysbyvectorizedexpressions. Thepresentprobleminvolvesa loopfor comput- ing the right-handside: for i in range(1, N): rhs[i] = (beta/dx**2)*(u[i+1] - 2*u[i] + u[i-1]) + g(x[i], t) This loopcanbe replacedbyavectorizedexpressionwith the followingreasoning. Wewant to set all the inner points at once: rhs[1:N-1] (this goes from index 1 up to, but not including, N). As the loop index i runs from1 to N-1, theu[i+1] termwill cover all the inner u values displaced one index to the right (compared to 1:N-1), i.e., u[2:N]. Similarly, u[i-1] corresponds to all inner u values dis- placedone index to the left:u[0:N-2]. Finally,u[i]has thesame indicesasrhs: u[1:N-1]. Thevectorized loopcan thereforebewritten in termsof slices: