Quantitive Methods in Hydro
Quantitive Methods in Hydro HYD 510
Popular in Course
verified elite notetaker
Popular in Environmental Science & Engineering
This 129 page Class Notes was uploaded by Dr. Wendell Schuppe on Thursday October 15, 2015. The Class Notes belongs to HYD 510 at New Mexico Institute of Mining and Technology taught by Staff in Fall. Since its upload, it has received 29 views. For similar materials see /class/223633/hyd-510-new-mexico-institute-of-mining-and-technology in Environmental Science & Engineering at New Mexico Institute of Mining and Technology.
Reviews for Quantitive Methods in Hydro
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/15/15
ProblemHSil 10f3 file CDocum ents20and20SettingsravindIaDesktophtm 1Prob1 Contents This is the backgroundsolution for a particle tracking problem Given a 2D velocity field uv in xy space and a starting location and time find the particle location at some later time Programmed by John L Wilson For Hyd 510 New Mexico Tech October 17 2008 DOES NOT CONTAIN THE PARTICLE TRACKERS The steady velocity field is u2cxy vca 2x 2y 2 for an incompressible irrotational fluid This is a quotpotential flow fieldquot ie satisfies certain conditions There is an exact solution for a streamline in this potential flow field It given by the parametric equation see our calculus review 28 psi ca 2 x xA3 3 yA2 x where psi is the quotvaluequot of the streamline Particle pathlines should follow a streamline consistent with this parametric equation Given a particle starting location this equation is solved for the value of psi along the streamline Solving this equation for y as a function of x and psi ysqrta02xA23 psicx This may be useful in plotting the streamline For this exercise the domain of interest is 0ltxlt5 0ltylt5 The particle tracking algorithem is based on the kinematic equations dxpdt u dypdt v where xp and yp are particle locations and t is time Two types of particle tracking are used a forward Euler method and a 4th order Runge Kutta scheme Definitions x and y are Cartesian coordinates m xp and yp represent particle locations m xp0 and yp0 represent starting locations m t tf and t0 represent time final time and starting time s uv are xy direction velocities ms Clear the screen and memomry c ear c c Hardwire velocity field parameters a2 cl Hardwire particle starting loction and time xp09 yp06 t00 Don39t yet constrain final time tf tfTnf 10192008 841 PM FmblanHSJ le c DucummtsmuandmusmngsmmnmDemmemPmm qume vebcmy e d and sheamhne xhnspscel yl yl h llnspsc l yl rlm xx meshgnmxm memy Wc a z ifzr lfz qulverlxy lyuyvh xllmHDJDNF yumlma n hold on Determlne and plot the screamllni passlng through 1 dd xhnspscel yl yl p5 a z xpEHLXpEI CHCPWPEI ZNXPD 2H A215lPs11gtlt1 Eu echhemeguesheve pmuxpwnw Undeilned fumclun or varlahle xp Error 111 gt hummus a an pmuxpwnw R unge Ku a Scheme gues have 2mm lUIDHUUHAIFM ProblemHSil file CDocum ents20and20SettingsravindIaDesktophtm 1Prob1 lotxpyp39x39gt add a legend NOTEs I found it best to use a while statement using it to iterate as long as the particle is within the desired domain hasn39t left yet the number of iterations is not too large and optional the time hasn39t exceeded some prescribed time useful if the time step is varied I also found it useful to write a function mfile subprogram that I could feed particle location xpyp and have it return the velocity to me This was particularly valuable to shorten the K scheme code D dpdpdpdpdpdpdpdpdpdpdpdpdpdpdpdp39o old off Published WW7 MATLAB 7 5 3 of3 10192008 841 PM New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Partial Differential Equations PDEs This is new material mainly presented by the notes supplemented by Chap 1 from Celia and Gray 1992 ito be posted on the web7 and Chapter 12 and related numerics in Chap 21 in Kreyszig Fundamentals of Partial Differential Equations We ll rst examine the motivation for studying PDEs then examine their nature and classi cation and nally talk about various solution methods Recall that the basic attribute of a PDE is that it has two or more independent variables In most applications these represent time and space We have one time dimension and up to three space dimensions These independent variables can be present in any combination but at least two of them must be present to be considered a PDE In presenting this material we usually represent time by the symbol t and space by conventional Cartesian coordinates x yz However as you see in the reading from Celia and Gray 1992 we can treat all four as independent variables without worrying about exactly what they represent Thus in a generic approach trepresents an independent variable which could be time space or even something else Since by convention the symbol y is now an independent variable we can no longer use it to represent the unknown dependent function In any one application we have a speci c symbol for the unknown such as p for uid pressure h for hydraulic head C for solute concentration T for temperature and so on In the mean time we need a generic symbol for the unknown to use when we are not examining a speci c application Celia and Gray and your text book both adopt the convention using it as the independent variable In general in threedimensional space and where t can represent time we have Independent variables x yz and t Dependent variable u uxyzt Of course to be a PDE there need only be two independent variables and that is where we start Consequently we examine problems where uxt or ux y As we examine things generically the symbols for y and t can be freely substituted for each other We are also examining PDEs that have both lst or 2quotd order terms or both in them as well as the possibility of zero order terms Most will have 2quotd order terms Motivation Most of the PDEs that we encounter represent a intensive conservation equation that is the conservation of mass energy or momentum per unit volume per unit time Each of these conservation or balance equations has three components which can be written out in simple terms for a hypothetical control volume as Time rate of change of storage Net ux rate in Net source rate 167 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology The units would depend on what is being conserved If mass then the rst and other terms have units of kg m39 s39 If heat energy the units would be J m393 s39l while if momentum they are kg m39 s39z Of course one could multiply through by a constant and get some other consistent units For example if the space domain is one dimensional we often multiply by a unit area m2 Or if preserving uid mass for an incompressible uid containing only a dilute concentration of solutes you might divide by uid density and end up with an equation preserving uid volume instead This latter approximation is common in hydrology The time rate of change term appears as something like amt or S auar where S is any storagelike coefficient For uids it represents compressibility effects for heat transport in a uid or solid it is heat capacity and for solute in a porous media it is porosity There are typically two kinds of uxes that concern us To explain this let s focus on conservation of solute mass but realizing that these concepts also apply to uid energy and momentum The rst ux is that due to solute advection with the owing uid solute concentration C times the uid velocity v or in other words ux VC This leads to a rst derivative contribution in the conservation equation Why The conservation equation is concerned with changes in concentration and the advective term speci cally deals with spatial changes caused by a owing uid Using a Taylor series approach we take the derivative to get at the rate of change eg 0c 6VC6x 1 That is why the 1st derivative ends up in the model The second type of ux is due to diffusion or something that has the same type of mathematical representation a socalled homology Molecular diffusive ux is modeled by F ick s Law in which the ux is proportional to the gradient of 39 The l l 39 quot f 39 D is called the molecular diffusion coef cient and in 1D the ux D 6C6x Again we use a Taylor Series approach to get the rate of change of diffusive ux eg 0c 6 D SC13x 6x D 32CEx2 if D is constant to insert in the conservation equation Many transport processes have this behavior such as uid momentum diffusion D kinematic viscosity Newton s law and thermal energy diffusion D thermal diffusivity F ourier s Law Other upscaled transport processes have a homologous behavior such as uid turbulent eddy diffusion D eddy viscosity porous media solute dispersion caused by spatially varying velocity D aquifer dispersion coef cient and aquifer ow D hydraulic diffusivity KSS where K is hydraulic conductivity and SS is specific storage Below keep you eye out for one or both types of uxes 7 advection and diffusion In our adapted notation the x direction uxes are given by vu and D wax which appear in conservation equations rate of change of ux as something like 6vu6x and D alt3x2 where v is velocity and D is the diffusion coef cient The net source term on the righthandside RHS of the conservation equation refers to sources sinks The sources and sinks can include rstorder terms that are proportional to state variable u and zeroorder prescribed terms not depending on u Some common generic PDE examples of relevance to hydrology are Bu flu D 2 6t 6x 0 The 1D diffusion equation uxt l 1 In our later development for many of the examples we assume v is a constant Thus 6vC8x vaCax 168 New Mexico Tech Hydrology Program Bu Bu v 0 The 1D advection equation uxt 6t 6x where v and D are parameters The 1D advectiondiffusion equation uxt Z Z a 1 a 1 The 2D LaPlace equation uxy 6x 6y 2 Z a 1 46 fxy The 2D Poisson equation uxy 6x 6y 2 Z Z a 1 a L a 1 The 3D LaPlace equation uxyz 6x By 62 Z Z Z 6 v6 v 61 1 6321 6 D 6 1 0 uxyt 6t 6 y 6y 6x W 6x6 W By A 2D advectiondiffusion equation 6u 6211 zu zu 0 The 3D heatconduction e uation u x zt at 6x2 Byz 622 q y Hyd 5 Quantitative Methods in Hydrology 2 3 4 5 6 7 8 where V vy Du ny and Dyv are parameters Notice in equation 7 we have a second order so called cross derivative term involving both x and y The presence of crossderivatives affects the choice of solution method Also notice that one of these equations has four independent variables two have three independent variables and the rest have two Each of these examples has been used to model solute movement and heat transfer for an appropriate conceptual model Equations 1 4 5 6 and 8 are also used to model groundwater owz Equations 3 and 7 are used to model flow in the vadose zone where v represents the in uence of gravity We could add rst order terms ax u and zeroorder forcings like that in the Poisson equation to expand the population of hydrologically relevant equations even more 2 You might wonder why we don t include the equation for uid flow in a river or the atmosphere The most general form of these equations is called the NavierStokes equation representing Newton s second law It is more complicated than the equations here and highly nonlinear Recall Newton s second law the rate of change of mom entum equals the sum of applied forces lts nearest relative above is the advectiondiffusion equation 3 The storage advection and diffusion terms of 3 would then represent the time and space rate of change of momentum You would add forces to the right side as net sources of momentum typically we add gravity and other body forces In any event our study of solutions of the advection diffusion equation provides the foundation for solving the NavierStokes equations and other more sophisticated models 169 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Notice that our list doesn t include a wave equation like 2 Z ozu 6 i 0 9 6 I 6x This looks similar to the LaPlace equation 4 but the difference in sign causes a signi cant difference in behavior modeled and in the nature and method of solution Your textbook focuses on this equation and it is also discussed in Celia and Gray 1992 It has interesting properties and wide application e g vibrations water waves ood waves seismicity etc It is just not so important or applied in hydrology as the other equations major exception using geophysical tools for characterization You can study this equation in your text Chapter 12 and in more advanced classes We skip it here Example Suppose a toxic solute is accidentally released to a river Shortly downstream it will have mixed transversely across the river due to turbulence Let C represent this crosssectional mean solute concentration The solute then advects downstream with the mean river velocity and disperses in the longitudinal direction along the river axis due to velocity differences across the channel velocity is lowest near the bottom and banks and highest toward the center of the rivers solute moves faster near the center causing an apparent spreading along the stream axis We model this as a socalled longitudinal dispersion process 0 10 Where v is the mean velocity and D is the longitudinal dispersion coefficient There are two independent variables x space and t time The highest order derivative in x which de nes it as a second order equation in space represents the longitudinal dispersion The first order derivative in x represents the advection while the first order derivative in trepresents the time rate of change of solute storage This is a commonly used first cut model of solute transport in a river3 3 Note that tracer tests in stream 5 designed to elucidate interchange between the stream and underlying sediments start with this model and then add a first order exchange term This is currently a hot topic l70 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Example 9 Water as a liquid moves through the vadose zone in response to gravity and gradients of pressure You may i I recall that the vadose zone has pore spaces lled with Z both air and liquid water The water pressure depends on the water saturation and related capillary forces Because the soil is only partially saturated the pressure is negative due to capillarity If the soil is uniform in composition capillary pressures are most negative where the soil is dry and most positive where it is wet A onedimensional vertical model through the vadose zone balances water storage measured by saturation ux due to gravity and ux to capillary diffusion the wicking of water from wet to dry areas As a PDE it can be written as 66 a 66 E ED6EK6 O 11 where z is positive down in this model and 9 is the socalled volumetric water content It represents the proportion of the space lled by water If S is the degree of saturation and n is porosity then 9 Sn D is the socalled soil moisture dyfusivity and K is the saturation dependent hydraulic conductivity As you can see they are nonlinear functions of 9 We will linearize l 1 here making D and K constants 0 12 As you can see this equation describing in ltration in the vadose zone is an advectiondiffusion equation The advection is due the gravity and the diffusion is due to capillary wicking The non linear version is simply a nonlinear advectiondiffusion equation The Nature and Classi cation of PDEs Notation To avoid vector calculus notation Celia and Gray p 9 1992 and your text introduce the commonly used and convenient subscript notation for partial derivatives Thus the generic PDE 2 6 V 6 D a 1 2 6t 6x 6x 0 13a becomes uvux7Duxx0 13b We use that notation here and in what follows 1st order equations Example Consider the 1D advective equation New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Bu Bu v 0 14a 6t 6x or uvux0 14b where x is space and tis time This is a common lst order PDE in hydrology usually encountered where u represents temperature or solute concentration Let s assume that the solution for this equation is in the spacetime domain 0 lt x lt L and 0 lt t lt T Then the solution to equation 14 describes how u varies in time and space within this domain but doesn t give us the actual values To get those values we need additional conditions As this PDE is lst order in space we need one boundary condition BC and as it is lSt order in time we need one initial condition IC Let s assume that v is velocity Then we need to prescribe one upstream BC If v is positive then we prescribe the BC at x0 The IC is prescribed at PO You can pick any other upstream location for the BC and any other time for the IC The in uence of the BC propagates into the solution at the speed given by v That is it propagates as a front 4 Those spatial locations where xgt Vt haven t felt the BC while those with x vt have The front propagates along a characteristic curve on the x t plane IT characteristic front I Hasn t felt the BC 00 In fact all portions of the solution propagate along characteristics5 What happens if the BC at x0 changes with time that is if u0t is not constant Then at each time tthe appropriate boundary value of u will propagate into the domain along its own characteristic 4 The term front suggests more than one spatial dimension Indeed there are 2D and 3D versions of 14 For example a 2D version is ut vX uX v u 0 The term front is commonly used in l 2 and 3D 5 One method of solution commonly used especially for this lst order problem is the method afcharacteris cs This concept is widely used even when solving with other methods l72 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Characteristic along which BC I T at t3 moves characteristics time changing BC t3 I2 Characteristic I along which IC 1 at xb moves 0 gt xa x L x b spatially dependent IC Problem 1 Suppose there is an initial spatial distribution of ux 0 the IC but that the BC is zero u0t 0 What would happen to the initial distribution What in uence would the BC have Below is a sketch of this situation where we ve kept some typical characteristics to help us Draw the distribution of u at two later times t1 and t3 given the IC shown on the sketch characteristics IT t3 t2 t1 0 0 xa xb L x Problem 2 Now suppose that the IC is ux 0 0 but that the BC is nonzero with the temporal distribution shown in the next sketch How would that BC propagate across the domain Draw it at two times when it has rst reached xa and then x6 l73 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology characteristics BC u0t gt 0 xa xb L x Page 11 0f Celia and Gray CG formally derive the idea of characteristics for this simple problem Let s review that here Characteristic Curves A characteristic curve is a curve along which a singularity will propagate De ne an arbitrary function gaxt continuous in independent variables x and t A constant value of o de nes a relationship between x and t representing a curve in x t space Assume that ux y is de ned along the curve of constant a Of course it is not itself constant along this curve Then the derivatives ux and u are de ned on the curve and are related by the derivative of it along the curve a t amt const Du ut D I uX Dx x 0 CG solve for the characteristics as ltoltox x0 vr r 15 where x0 and t0 are arbitrary constants think of them as starting points or IVs Note that go gpx Vt x0 vt0 changes only when x vt changes There is a different characteristic for each value ofx vt Thus the family of characteristics is de ned by x Vt constant 16a or 174 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology x x0 vt t 0 16b These are the curves lines we plotted in the gures on the previous page 16b is the referred to as the characteristic equation for 15 because it gives a set of characteristics curves for different values of the constants6 Complete knowledge about it along one of the characteristic curves in NOT suf cient to know about it along any other curve Suppose instead we had complete information along another curve that intersects crosses the characteristics Then we would have a value of u at one point along each characteristic and therefore everywhere along all characteristics and throughout the domain We illustrated this by looking at the propagation of initial and boundary conditions each lling half the domain with information We need to know both to ll the entire domain 2quot 1 order equations Celia and Gray p 15 1992 write the general 2quotd order PDE for uxt Au zBu Cu Du Eu Ru x t 177 where x and t are independent variables I could be space as well as time just replace the symbol t by y Note the crossderivative on the left side Also note that the forcing on the right side which now includes the 03911 order 2 u terms The behavior of the solution and the choice of solution method depends on the presence sign and value of the coef cients A B C D and E This suggests a classi cation scheme for 17 based on these coef cients and the behavior of the solution Your text presents a version of this classi cation on p 551 but the version in Celia and Gray CG is both more complete and easier to grasp We won t bother here with the derivation see CG pp 1419 which employs Cramer s rule but rather talk about the result and then examine the implications which are signi cant There are three classes of equation 14 based on the values ofA B and C the coef cients for the second derivative terms These classes are p 18 CG lfBziAC gt 0 the equation is hyperbolic lfBziAC 0 the equation is parabolic 18 lfBziAC lt 0 the equation is elliptic 6 Although there are an infinite number of characteristic curves each conforming to a particular set of constants 15 is said to have one characteristic 7 Note the D in this equation represents any coefficient of the lst derivative wrt I Don t confuse it with the diffusion coefficient which is associated with the second derivative wrt space l75 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology The classes are based on solving for characteristic curves similar to the discussion for the ISI order equation Let s illustrate by examining individual cases using domain boxes that diagram the behavior This will also bring in boundary and initial conditions When the PDE in 17 is 2quotd order in x or t we will need two boundary conditions BCs or initial conditions le but when it is only 1st order in twe need only one initial condition IC Finally keep in mind that in hydrologic systems we seldom have problems in which B is nonzeros Elliptic Problems The LaPlace equation 4 is typical of elliptic equations 2 2 a i a 19 6x 6y Comparing to the general equation in 17 where t is replaced by y we see thatA and C are non zero and equal 1 while B D andE are zero and so is R ConsequentlyB iAC l lt 0 and the equation is classi ed as elliptic Typical applications are9 heat conduction uT temperature 0C groundwater ow uh hydraulic head m and molecular diffusion uC soluted concentration kgm3 Let s focus on heat conduction In heat conduction equation 19 is used for example to model heat transfer in a 2D plate or one of its geophysical analogies10 Key Cool Warm T06 V x Panoramlc v1ew Top v1ew Heat conduction in a 2D plate 8 Wave equations are typical of cases in whichB is nonzero They have eigen solutions with multiple modes characteristic frequencies that we will ignore this semester Your textbook concentrates on these problems in Chapter 9 As quasisteadystate processes see next page for explanation 10 lnsomuch as groundwater ow and heat conduction are homologies we borrow heat conduction solutions to represent groundwater flow Carslaw and Jaeger s 1959 Heat Conduction in Solids is one of the most cited publications in the groundwater literature This literature is also the source of the solution used by CV Theis to represent transient groundwater drawdown around a pumping well l76 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology We need four boundary conditions two in each direction as the problem is 2quotd order in both x and y These are labeled BCs 14 in the sketch We can use 151 2quotd or 3rd type BCs We introduced these on p 122 of these notes In this context a general boundary condition is written as k1 u kzun K 20 along a boundary where the k s are parameters weighing the relative contribution of u and its gradient uquot at each boundary The gradient uquot is the derivative taken normal to the boundary un ux for BCs 2 and 4 in the sketch and uquot uy for BCs 1 and 3 The K s are values Along first type or Dirichlet boundaries k1 1 k2 0 and K is the known numerical value u on the boundary Along second type or Neumann boundaries k1 0 k2 1 and K is the value of uquot on the boundary Finally if both k1 and k are nonzero we have a third type or Cauchy boundary also called aRobin boundary as in CG All three boundary condition types are used to solve equations of the type given in 19 and in 1 and 3 8 on p 168 ofthese notes In our heat conduction example if all four boundary conditions are known temperature T 1St type boundaries then the problem is called aDirichlet problem In the sketch the temperature is indicated by the shade of gray darker gray being a warmer temperature Can you infer the likely boundary types and values from the temperature pattern The characteristics for an elliptic problem with B2 iAC lt 0 are real In short elliptic problems instantaneously and completely feel a change in boundary condition or a sourcesink which we could add to the right side of 19 In this example x andy represent space Suppose a third independent variable trepresents time Now further suppose that one of the boundary values say a 1St type boundary changes in t In other words a prescribed boundary temperature value on one of the boundaries changes Let s make it a step change as shown in the sketch at left Before tt1 the prescribed temperature on the boundary is To while afterwards is T1 T T To T1 gt t I Step change left and continuous change right in boundary temperature Then With the step change on the boundary the solution in the interior immediately adjusts to the new boundary value There is no delay There are in essence two solutions one before the change and one after Let s make it more interesting Suppose that along this Dirichlet boundary Tt varies continuously in time as shown in the sketch at right then the solution instantaneously and completely responds to this changing boundary It s as if the solution was composed of a series of continuously changing steady states that are always in equilibrium with the boundary conditions Let s visualize this with using a domain box that diagrams the behavior 47 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology A y No real characteristics in the elliptic problem Feels the BCs immediately 0 gt 0 L x Parabolic Problems The diffusion equation 1 is typical of parabolic equations Bu flu 0 21 6t 6x2 where x is space and t is time Comparing to the general equation in 17 which has both x and t we see that coefficientsl and D11 are nonzero A D D 1 while B C andE are zero and so is R Consequently B2 iAC 0 and the equation is classi ed as parabolic But the nonzero value of D implies a secondary characteristic see p 18 CG that is interesting Typical applications are similar to the LaPlace equation except now they involve time transient behavior heat diffusion uT temperature 0C groundwater flow uh hydraulic head m and molecular diffusion uC soluted concentration kgm3 Again we focus on heat conduction In heat conduction equation 21 is used for example to model transient heat transfer in a 1D bar or rod or one of its geophysical analogies Key I4 I Cool t3 I3 Warm t t2 I3 BCl BC2 Q n TOW t00 Transient heat conduction in a 1D rod 11 Don t confuse the two differentD s 12 A typical analog is vertical heat conduction through soil or though an aquitard Homologies include vertical water leakage of an aquitard with storage in the aquitard horizontal capillary diffusion of moisture in an unsaturated column l78 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology We need two boundary conditions one on each end of the rod as the problem is 2quotd order in x These are labeled BCs 1 and 2 in the sketch We can use 151 2quotd or 3rd type BCs We need one initial condition Tx0 as the problem is lSt order in time I There is one real characteristic for a parabolic problem with B2 iAC 0 The characteristic is pp 22 CG go 01 to or since g9 constant along a characteristic pro 0 22 The sketch below shows the characteristics for this parabolic problem The characteristics are horizontal lines parallel to the x axis Disturbances at the ends of the rod propagate instantaneously to the interior However experience suggests that this doesn t actually happen The primary characteristics don t tell the full story One has to look at the secondary characteristics associated with the 1st order term t to 0 time changing t3 t2 t1 0 gt 0 L x The secondary characteristics have the form x 7 x0 0 23 These are shown as dashed lines in the next sketch or transient diffusion of a solute across a barrier which can range in scale from an intrapore bacterial polysaccride to hundreds of meters of shale l79 New Mexico Tech H d 510 Hydrology Program Quantitative Methods in Hydrology Characteristics x x00 spatially dependent IC In short parabolic problems remain sensitive to initial conditions which propagate along the dashed character1st1cs and also instantaneously responds to a change in boundary condition however slight that response might be which propagate along the solid characteristics From your CG For point P in Figure 18 the domain of dependence is the region 0 lt CL lt L below P the shaded region Information propagates toward P from both ends of the rod along a horizontal characteristic Because the characteristic is horizontal information from both ends of the rod is predicted to propagate instantaneously to any point on the interior of the rod although the instantaneous effect at points far from the boundary would be small However the temperature at a point such as P is also in uenced by the previous temperature pro les in the rod Thus the secondary characteristics which emanate from t O are also signi cant in calculation of the temperature Because only one type of characteristic emanates from the z axis into time only one initiai condition must be speci ed This condition must involve derivatives of order less than the highestorder time derivative appearing in the differential equation Thus the appropriate initial condition is speci cation of T all along the rod The primary characteristics are horizontal Therefore information may propagate into the rod from both ends and boundary conditions must be speci ed at both ends Because T is speci ed as an initial condition the boundary conditions may be either rst type second type or third type The suf cient boundary conditions are presented in Figure 19 The parabolic problem studied here is a mixed initial in time and boundary in space value problem i x 2 0 x L Figure 18 Domain of dependence for solution of equation 0342 at point P Solid lines are primary characteristics from equation 1346 and dashed lines are secondary characteristics from equation 1347 l80 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology I Specify Specify T T or g T T 9 91 QT 6T HT 5T 1 z 0 x 0 x L Specify T Figure 19 Boundary and initial conditions for parabolic equation 1342 solved in region 0 lt a lt L and t gt 0 The Advection Diffusion Equation Consider the threedimensional advectiondiffusion equation with rst order decay CG l3639 see 7 on p 169 of these notes for a 2D version but with directionally dependent diffusion and Without the lSt order decay term This extends us in two ways First it adds rst derivatives in space the advection to the parabolic equation in 21 and it takes us to four independent variables three in space and one in time This equation is typical of solute transport but also models heat transport due to advection with owing uid along with heat conduction For a solute the decay could be due to radioactive decay or biologically mediated transformations The equation is 1396 in CG39 modi ed for our notation du du du du 6214 6214 6214 ku 0 23 6x2 322 622 v V vz at 6x y 6y 62 Where u is concentration the v s are components of velocity D is the diffusion coef cient and k is the rate constant This equation is parabolic and has primary characteristics 0 t7 t 0 24 but it also has very important secondary characteristics p 39 CG v v vt t vx x0 vyy y0 vz 20 0 25 If the v s are not constant then these characteristics have curvature in 3D space If they are constant they describe hyperplanes analogous to the 1st order characteristics discussed on p 173 these notes and shown in 16b ie vt t0 vx x x0 0 For example if we write 23 for one spatial dimension we have 2 a va Da ku0 26 at 6x 6x with characteristics illustrated by this sketch l8l New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology time changing BC Characteristics vxt t 0 xxo0 gt 0 L x 182 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Numerics for FirstOrder ODEs Background and motivation JLW notes here also read the text s 211 and 212 An additional reference available from the TA is s 8187 iespecially s 8185 of Gilat and Subramaniam Numerical Methods for Engineers and Scientists Wiley 2008 Analytical solutions to differential equations are sometimes hard to find The effort can be worth it if the solution is easy If the solution is complicated it can still be worth it if the solution will be used over and over again or if the symbolic form of the solution is selfrevealing and perhaps could be published On the other hand if the solution is to be used once or only a few times if only numerical values of the solution are sought or if the analytical solution is too complex or even impossible a numerical approach is appropriate Consider the IV problem y39fxy yx0yo 1 where x0 and yo are given and the solution is required and exists on some nitesolution closed interval ab Depending on the numerical method used to solve 1 the solution is obtained at a nite number ofpoints xquot n 0 1 2 3 N in the closed interval ab For an initial value problem x0 61 represents the location of the IC and xNb represents a stopping point If the solution is asymptotic to some limit b is simply selected to be large enough such that y closely approaches that limit The solution consists of values yquot for n 0 1 2 3 obtained at each of the points xquot where the rst value yo is the prescribed IC The solution then consists of the N1 Pairs x0 yo 961 M 962 yz x3 y3 xN yN Figure 1 Numerical discretization of an ODE Looking ahead for motivation we ll see that a system of ODEs consists of solving the simultaneous system of equations Y39 fxY YOCO yo 2 where y is a column vector of unknowns say y 141 yb yo yd where the subscripts a b c represent the different unknowns at each scalar point x and the prime indicates a vector transform New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology to turn the row vector into a column vector We use the same methods to solve vector oale 2 that we use to solve scalar ode 1 The equations in vector equation 2 may represent a number of different strongly coupled processes where the value of one unknown depends on the value of another An example is strongly coupled uid heat and mass transfer where the heat transfer coefficient depends on the ow rate and the ow rate depends on the heat transfer say through uid viscosity or density The other application of 2 is in the numerics that are used to solve PDEs In this case the independent variable x represents time ie x t and the unknowns in vector y represent spatially distributed values of a state variable eg solute concentration taken at a nite number of discretized points node points in space An example is a groundwater ow model such as that programmed into the applied code MODFLOW where in 2 y is a vector of nodal values representing spatially distributed hydraulic heads The nodal solutions are calculated at times t1 t2 t3 t3 by solving 2 with time t replacing x and hydraulic head h replacing y Approaches to initial value problems JLW There are a variety of issues to consider when setting up and interpreting a numerical solution to 1 How big do you make the step size Ax between the points xn1 and xquot Given a value of the unknown yquot at point xquot how do you estimate the value yn1 at the next point an What is the error in this estimation Is the estimation numerically stable Taylor Series and numerical approximation Most of the methods we ll consider are based on the Taylor Series or something similar We use the notation of Figure l The Taylor Series for a function y at location xn1 given the value and derivatives of y at location xquot is sz Ax3 yn1ynAxy39xn 2 y xn 3 y 39xn 3 where Ax xn1 xquot This is exact If however we neglect the higher order terms we get the secondorder approximation M H 2 y 5 4 The last term on the RHS respresents the secondorder local ym er Axy xn where xquot Sinn Axx truncation error see below for this approximation which we encounter by dropping the higher order terms n1 39 We are interested in solving the ODE y39 f x y thus we seek a numerical approximation to the derivative on the LHS of the ODE The Taylor Series provides such an approximation if we solve 4 for the first derivative ym yn Ax 39 x slo e m H 5 y n p 2 y 5 We now insert this approximation into the ODE to get the nite ah erence model ym yn Ax H m x n 6 2 y 5 f n y where we take the RHS of f at the same point xquot around which we write the Taylor Series Solving for yn1 gives New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology sz ymey m xmyn 2 yquot 7 which we usually write as yin eyi AxfxiynOAx2 8 to indicate the order of the approximation This particular approximation is called a forward or explicit Euler approximations see pp 886888 oftext Note that the explicit method handles nonlinearity by explicitly evaluating f at the previous location xquot and previous value of the unknown yquot yx exact solution local truncation error numerical solution Figure 2 Euler s explicit method From this example we can extrapolate to a general form to represent other approximations ynl mFxn xnl yn yn1OAxk 9 where k is the order of the approximation The right side function F can include a contribution from the current unknown ynH It also implicitly includes some combination of 2 and m For example in the explicit Euler method 8 there is a potential contribution from xquot xn1 yquot and 2 to the right side but there is no possible contribution from yn1 or m For a given ODE different numerical methods yield different functions F and different approximation orders k Errors The numerical methods used to make the estimation are subject to both roundoff errors and truncation errors Roundoff error is due to the limited word length of a computer while truncation error is due to the numerical approximation nature of the method used to calculate the solution as in 9 That is roundoff error is due to the computer while truncation error is within your numerical algorithm You can select better numerical methods see below that reduce truncation error increase the value of k in 9 but only at the expense of additional computational effort Or you can reduce the step size Ax to reduce truncation error and keep a 5 Both of these terms are commonly used New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology simpler method However there is a limit to the effectiveness of reducing step size However when the step gets too small roundo error dominates and the solution is worthless Global V local truncation error Truncation enror consists of two parts Local truncation error is due to the application of the numerical method to a single step See 7 above to observe that the explicit Euler method has a local error of order sz when measured as to how well the unknown y is approximated Propagated or accumulated truncation error is due to the accumulation of error from previous steps Together they produce global total truncation error in the solution The explicit Euler method has a global error of order Ax the number of steps is proportional to 1 Ax which you multiply by the local error sz in each step The order of approximation for global error is almost always lower order than the local error Balance error The error measures above look at how well the function y is approximated But we are hydrologists and we are particularly interested in how well the balance of mass momentum or energy is represented by a numerical approximation Consider again the ODE y39 fx y Rewrite is as a balance equation y fx y 0 Ifwe replace the true solution by a numerical approximation then it can t possibly balance at least not perfectly We can only hope to get close Rewrite the balance as y39 fx y 8 10 If we have the exact solution for y then the outofbalance error 8 will be zero But when we replace y by its numerical approximation 8 will be nonzero How much different from zero will it be For the explicit Euler method you can infer from 6 above that its local balance error is only of order Ax and that the global balance error is even worse Aside when solving nonlinear versions of y fx y 0 especially when involved with simultaneous solutions as in 2 we often resort to Newton s method p 790 of text searching for the value of y that takes 3 in 10 close to zero That is we search for the root of y fx y 0 Numerical stability Many of the numerical methods used to solve l are unconditionally stable Some methods we like to teach which are used by beginners and which lend themselves to parallel computation are conditionally stable an example is Euler s forward or explicit scheme presented above Make the Ax too large and they blow up Of course some ODEs are themselves unstable we are referring here only to instabilities introduced by the numerical method There are also intrinsically unstable numerical methods for a given ODE that we will avoid Step size If uniform step sizes Ax are used mathematicians usually refer to the step with symbol h or hAx constant In many applications we allow the step size to vary either prescriptively or adaptively In prescriptive declarations of step sizes the user decides where ner resolution is needed In adaptive methods an algorithm is supplied to make this decision as part of the solution In any event the step must be small enough to minimize truncation error to an acceptable value but not so small as to encounter unacceptable roundoff see above Prescription of variable step size In many cases you can anticipate the general functional nature of the solution that is you know something about whenwhere y is going to vary rapidly with x where the action is and whenwhere it is not In this case you can prescribe a variable New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Axnxn1 xquot increasing the density of the points xquot where the action is in order to reduce truncation error there and reducing it elsewhere Two examples are groundwater flow and heat transport both of which respond exponentially to change When a new forcing is added you can anticipate that you ll need smaller Ax s after the change but can enlarge them latter For example many groundwater codes prescriptively reduce the time step after a well is turned on and then let the time step increase logarithmically afterwards Adaptive step size In other cases you don t know where the smaller time steps are needed after all they are a function of the asyet uncalculated solution Suppose that you want the numerical method to figure this out and automatically reduce Ax wherewhen the action is and then increase it away from the action A variety of adaptation schemes are available See p 889 of text for a general overview which is illustrated with the explicit Euler method These adaptation schemes are usually based on nding the rate of change of the ODE y39 f x y with respect to wrt the independent variable x Using the chain rule y f39 fyy39 fyf 11 They then select the next step size Axnxn1 xquot such that the truncation error does not exceed a certain tolerance Explicit V implicit methods The difference between explicit and implicit methods is in how the solution is calculated at each time step Explicit methods use an explicit orforward formula for calculating yn1 at point xn1 given yquot and at point xquot In an explicit formula the right side of the numerical equation has only known quantities That is the ODE y39 f x y is replaced through the general numerical approximation in 9 by the particular approximation ym Fxnaxwyn0mk 12 The formula is explicit in that the variables in the right side function xn xn1 yquot are all known The unknown yn does not appear on the RHS For this reason explicit methods handle non linearities by simply evaluating their contributions to F using information already available at xquot The explicit Euler approximation in 8 is the simplist example of this With implicit methods the unknown yn1 appears on both sides of the equation just as in the general model 9 ynl mFxn xnl yn yn1OAxk 9 For example suppose that in the Euler method we took the Taylor Series around xn1 instead of around xquot Then 4 would become ml 2 y e ym Ax y39xn1 y 13 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology where again x S 5 S x Ax cm1 and we still have Ax xn xquot Note the change in sign on the rst derivative odd term since we are now going backwards The nite difference model becomes yn1yn u N 7m 2 y fxn1yn1 14 Herer we take the RHS value of f at the same point xn1 around which we write the Taylor Series And that point is now the same as the point xn1 at which we are solving for yn1 0r rm e y Ax fxi1yi1OAxZ 15 This is an implicit model because yn is now on both sides ofthe equation We call 15 a backward or fully implicit Euler method see p 896 in text Notice its truncation error is similar to the explicit scheme Also note that because yn1 is also on the right sides nonlinearities must be handled implicitly Aside most lV problems in hydrology have time as the independent variable Both implicit and explicit methods are used to solve these problems Most spatially distributed problems are boundary value problems BVP involving second order ODEs in space BVPs problems must be solved implicitly Single Step Methods The Euler schemes are examples ofa singlestep scheme Most hydrologic models use single step schemes so it is a good idea to master them In a single step scheme the numerical solution yn1 at a point xn1 is calculated from the known solution yquot at the previous point xquot The general algorithm to solve l is xnl m xn Axquot 16ab ym e yi An slope where Axquot xn1 xquot is the step size and slope is an estimate of the value of dydx in the interval xquot to xn1 The solution starts at the point x0 where the solution is known and marches forward recursively step by step n0 l 2 N see Figures 1 and 4 A yx y exact 7 solution L y Slope numerical solution yr quot i xn 4 p xnl x Ax Figure 4 Singlestep explicit methods New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology The difference between singlestep methods is in the value used for slope Methods become more sophisticated by improving the calculation of slope for example by subdividing the interval between xquot to xn into smaller pieces The order k of the approximation gets higher better with increased sophistication but the computation becomes more expensive because there are more arithmetic operations In explicit schemes the slope depends only on known values yquot and fquot while in implicit schemes it also depends on unknown values yn and m The explicit or forward Euler scheme 8 is the simplest singlestep method but there are other Euler schemes that are commonly used These include the backward scheme 15 and a mid point Euler scheme not in your text with theoretically less truncation error In any event Euler schemes are commonly encountered in most heat transfer solute transport groundwater and vadose zone codes The next most sophisticated singlestep class of methods are called Runge Kutta methods RK methods In these methods the slope is calculated from a weighted average of estimates of the slope at several points along the interval xquot to xn1 Higher order RK methods use more points MATLAB has Mfiles that use RK methods to solve ODEs and has several orders of RK methods available to choose from RK methods are typically used in pairs By comparing the solution at one order to that at another lower order you can estimate error and if needed decrease Ax RK methods are commonly used in hydrology for particle tracking in which the pathlines and travelresidence times in various hydrologic systems are calculated Multistep Methods In multstep methods the solution yn at a point xn is calculated from the known solutions at several previous points e g at three points yquot at xquot AND y at xn1 AND ymz at xnz The idea is that by using more points you can get a better idea of the trend of the solution These methods have not been widely applied in hydrology in the past but that is changing Also they are becoming common in sophisticated commericially available multiphysics and domain speci c codes like COMSOL Multiphysics and the CFD computational uid dynamics code Fluent which are finding more frequent use in hydrology Multistep methods are introduced in Section 212 ofthe text Predictor Corrector Methods Predictorcorrector methods use two formulas to solve an ODE at each step The predictor uses an explicit formula that is used to get an estimate of the solution yn1 Since this is an explicit calculation it only uses the known solution yquot at the previous step single step methods or several previous steps multistep methods Once the estimator yim is found a correction is applied Stiff ODEs Error terms involve the higher derivatives What happens to error when Ax increases If the error grows faster than the solution we have a problem the solution is unstable and can become nonsensical Such a system for which Ax must be restricted to smaller values and the physical systems they model are called sti Some solution methods can handle stiff problems like the backward Euler method others like the forward Euler method can only handle them by restricting Ax Later we ll examine multiple simultaneous odes Let s say the independent variable is time and the equations are linear for simplicity We ll discover another version of stiffness when the various odes have considerably different time constants Such stiff problems are harder to solve New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Euler s Methods Used extensively to solve a large number of hydrology problems including heat transfer uid ow in channels and overland groundwater ow vadose zone ow solute mass diffusion and advection dispersion etc Euler s Forward or Explicit Method See p 8990 above for derivation and pp 887889 in text for background and example Works on linear or nonlinear problems but is only conditionally stable ym e yquot N fwd0 0Ax2 8 exact yx solution yxn1 ym without additional local truncation error slope fxmyn numerical solutions ynH Figure 2 Euler s explicit method numerical and true solution Example from text Apply the Euler method to the following initial value problem choosing ll 02 and colnpulingyl 4 5 4l y 39 n W n Solution Here f39 y X 339 hence f39n ll39nl xquot i39n and we see that 3 becomes quotn1 7 02w w t 11 qu Mnhad Applied u 4 m Exlmplm I and Error i quotluv 7 U Fmquot N noun on quotmm mm Where an y y Euler s Backward or Fully Implicit Method See pp 9293 above for derivation and p 896 of text for discussion and example The method is unconditionally stable New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology ymi N y Ax fxmigtymi 0M 15 Evaluating the righthand side is no problem for a linear problem but can require something like a xedpoint or Newton iteration for a nonlinear problem The method is also particularly useful for stiff problems text p 896 Euler s Midpoint or the CrankNicholson Method This method is not mentioned in the text We review it below as a special case of the socalled general Euler 6method The midpoint method is higher order second order and stable and o en used in hydrology Stability Issues We ll review these in discssing the general Euler 6method below Modi ed Euler s Method a PredictorCorrector Approach This is called an improved Euler method or Huen s method in the text pp 889891 which includes pseudocode see below and an example Note that the rst equation is the prediction step and the second is the corrector step The prediction step is a simple forward explict prediction ym yo Ax fawn 8 ymi yl fltxlyigtltxmiyw As one would expect this has higher order local and global error see p 891 for proof The local error is of order Ax3 and the global error is of order sz In other words this method is a second order method Also note how this method handles nonlinearities The pseudocode is not 21 2 improved Euler Method Heun39s Method ALGORITHM EULHR l r quot0 ll Ni 39llii I iilllnmlupuiuxlllcxululimiullhciilllilxluluuplohluililquot 7 39 liiolliii ll l ml 4x II 1 i v r 1v it llmi illi pmiilom I illeilE mluliuii on the Hi ul Ii PM A lNl U39l39 lliiliiil lllucs i i wlcp size Ii Inlnllwr oi lk39ps iv UUTPU r Appminuiliun rm in tile mluliull l39H mll H i rquot H im VINICH U N l i din A lulu u Ail ii k2 slop Eull EULER New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Exact solutions to an ODE In order to better understand truncation error in the Euler method we need to examine exact solutions written in a similar way and compare to the various Euler approximations In this section we look at exact solutions In the next we reexamine approximate solutions and truncation error Keep in mind that in this section everything is exact even though we are representing solutions by in nite series In the next section we truncate those series to get the approximate solutions Consider a rst order ODE in the interval x0 S x 5xN It can be linear or nonlinear The equation is y x 000 1 where y is the unknown x is the independent variable f is a lnction of x and y y x is the derivative of y at location x This equation applies anywhere on domain along the xaxis where f is de ned Now consider 7 0 l 2 3 N7 l Nspeci c points xquot on the xaxis each separated from its n1St neighbor to the right by some distance Axquot xquot 1 7 xquot Note that unlike our discussion above we are now allowing the point spacing to vary Equation 1 applies everywhere in the intervale S x 5 xN and therefore it applies at each point xquot and its neighbor to the right xquot 1 Assume that 1 has exact solutions at these two points respectively yxn and yxn1 Where 1 is written for each point as y xn xmyxn lm 2 y xn1 xn1yxn1 xnl 3 On the right I ve introduced shorthand notation for the function f to indicate its x location not its dependence It is understood that it can depend on the value of both x and y at that location as explicitly shown in 1 Forward Taylor Series Expansion The solution at xn1 can be written as a Taylor series expansion around the solution at xquot As long as all terms in the expansion are retained the expansion is exact y xn yo m x x m xquot13xquot y39 xn 4 This is called a forward expansion It projects forward in x from xquot to xn1 Recalling the de nition of the increment or spacing Axquot this forward expansion becomes New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Ax 2 Ax 3 yxn1yxnAxn y xn 2quot y xn 339quot y x 5 We can now solve 5 for y xn which we ll eventually substitute into 2 yx we Axf y xn y 39xn 6 n 2 3 This expansion for the derivative at xquot is also called a forward expansion Once again it is exact When 6 is substituted into 2 we get y39xn yxquot1yxquot y x Axquot y 39x foe 7a 2 3 n 01 l l MHZ l l l n 2 y 960 3 y xn 7b which is a new but still exact version of 2 written as a forward difference and series expansion around xquot Expansions 5 and 7 will later become building blocks for nite difference numerical approximations when the higher order terms are dropped mix yx M Ax Since by 2 y xn 09 these expansions can also be written in terms of derivatives of the function f For example 5 becomes yo yxAx fx A f39x A fquotx 8 where e g f is the partial derivative of f wrt x And 6 becomes yxn1 yxn Axquot an u yxnTn Tfxn 3 f xn 9 When 9 is substituted into 2 we get a new version of 7 yxn1 WC AX Mil n T Tfmf 3 f xn fxn 10a or yxquot1 yxquot fx f39x Ag fquotx 10b The higher order terms of these series are nonzero only if the respective derivative wrt x of f is nonzero For example if f is a constant then only the rst term on the RHS of 10b is nonzero If f is a linear function of x and therefore the ODE 2 is linear then only the rst two terms are nonzero Etc Backward Taylor Series Expansion The solution at xquot can be written as a backward Taylor series expansion around the solution at xnl New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology NC yxn1x xn1y xn1 xnxn12 u 2 y xn1 3 3 11 xquot xquot1y39 x Substituting step size Axquot xn1 7 xquot this backward expansion becomes 2 3 Ax yxnyxn1Axny39xn1 y xm 339quot y 39xn1i 12 which resembles 5 but with alternating signs since x is going backwards and with the derivatives of y taken at a different point We can rewrite in terms of derivatives of the function f n Ax 2 M m3 yxnyxn1 Axn fxn1 2quot f39xm 339quot f xn1i 13 and then rewrite it for y xn 1 M m3 yxn1yxnAxn fxn1 2quot f39xn1 339quot f xn1 14 We can also get a backward expansion for the ode solving 12 for y xn1 which we then substitute into 3 The backward expansion is yx1yxn H M yxn1 Ax 2 y xn1 3 y xn1 and the ode becomes yx yx Ax Axf Ax fxn1 2 y xn1 3 y xn1 or yx yx Ax Axf Ax fxn1 2 fxm 3 f xm 17 which has been written in terms of derivatives of the function f These expansions at xn1 for the derivative and the backward difference are also called backward expansions Euler s Forward and Backward Methods We now develop approximate models for y and y at xquot and xn1 by dropping the higher order terms in the Taylor Series expansions We approximate 2 by a truncated forward expansion and 3 by a truncated backward expansion Notationally where yxn1 is the exact solution for y at an we write yn1 to represent the approximate numerical solution at the same point That is subscript notation is used to denote a discrete approximation Note that we hope that yn Z yxn1 but how well we meet this goals depends on the truncation error and the stability of the method Euler s forward method The exact forward expansion at xquot for the ODE 2 is given in 7b It can be rewritten as New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology yX1 WC f AX H x y 5 18 Ax 2 where the last term represents the rest of the series and 5 is somewhere in the interval xquot to xn1 That is this one term is exactly equivalent to the rest of the series if you know how to pick 4 If we neglect higher order terms then the no longer have an exact equation but a numerical approximation We write this as yn1 yn N Ax fn OAxn 19 By truncating the higher order terms in 18 the solution to 19 is not exact solution of 2 The second term on the RHS represents that truncation error Note the important notational change Where yxn1 is the exact solution for y at xn1y and satis es 18 exactly the variable yn1 is the approximate numerical solution at the same point and 19 is only an approximation The second term on the RHS of 19 represents that truncation error and is used to indicate that the order of local ODE error is Axquot That is it depends on the step size taken to the rst power For example if the ODE is one of mass balance then this represents the local mass balance error Often we don t write the truncation error term instead just writing the difference equation 20 The forward expansion for y from xquot to xn1 is given by 5 and 8 We approximate it by truncating the higher order terms to yield an approximate numerical solution for yxn1 Ax 2 y s y Ax fx 2quot yquotx 21 01 Ax 2 y s y Ax fx 2quot f x 22 We write this as y s y Ax fx0AX2 23 to indicate that the order of local truncation error for yn1 is Axnz and then without the truncation error just to indicate the solution to the forward difference 20 24gt The function f in the second term on the RHS is evaluated the old location xquot Euler s backward method We use equations 12l7 truncating higher order terms to nd the differences used in Euler s backward difference method From 16 or 17 we get the backward difference approximation to ODE 3 x 0Ax 25 or writing it without the truncation error 100 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology 26 The backward expansion for y from xn1 to xquot is given by 12 and 13 Approximating 13 we have yn E ym Axn fxn1 0Axn2 27 However we are interested in nding yn1 given yquot The expansion for this is 14 which we approximate as y s y Ax fxi10Axf 28 We can also get this from 27 It is usually written without the order of local truncation error ym E y Axquot f 9 29 Comparing the backward and forward methods we see that they have the same order of truncation error but that the function f in the second term on the RHS is different For the backward method it is evaluated at the new location xn1 Later we ll see that the choice of a backword or forward method turns on stability and computational effort General Euler 0 Method We can approximate the derivative y at any point xng OS 9 51 within the interval xquot to xquot 1 by a linear combination of the derivatives at the end points y39xigE19y39xii9y39xi1 30 where the derivatives on the right side are exact and are given by 2 and 3 i9 is a weighting coefficient We later replace the derivatives with their numerical approximations We also have the following special values of i9 and their associated numerical approximation 0 forward or explicit Euler Method l backward or fully implicit Euler Method 05 midpoint Euler Method 0 05 also called CrankNicholson when used to solve time dependent PDEs CDCDCD We ve already seen the results for 9 0 and 1 For 0lt i9 51 the method is implicit For 9 1 it is fully implicit Your textbook refers to the fully implicit Euler method simply as the implicit method It is the only implicit Euler method discussed in there The general method and in particular the midpoint method are not discussed However they are commonly used to solve the time derivative aspect of transient PDEs in hydrology so we discuss them here The exact ODE l at point xng is given by New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology y39xnafxng 31 We can write the LHS in terms of the weighted derivatives at the end points given in 30 We can also write the RHS in terms of f weighted at the end points fxng51 i9fxn 9fxn1 32 Note if we substitute 30 and 32 into 31 we get 1 9y39xn 9 y39xn11 i9fxn 9 fxn1 which satis es 1 exactly if we ignore the fact that 30 and 32 individually are not exact To proceed we need the numerical approximations to the derivatives in 30 and the values of f in 32 For example xquotfxnyxn where yxn is approximated by yquot in the numerical solution Truncate the series in 6 to develop a forward difference to approximate y Jen and truncate the series in 15 to get a backward difference to approximate y xn or v N yn1 yn v N ynl yn x and x 33ab y n n n1 n Both differences are of order Axquot and while the derivatives are at different locations x their differences RHS of 33 are the same Combine 3033 to get jf za m 6 u ea where represents the numerical approximation fquot to fxmyxn etc The numerical solution yn1approximating yxn1 is given by nuznAnb 9 n w As you can see this becomes the forward difference solution 24 when 9 0 and the backward difference solution 29 when 9 1 Its easy to shown that its local truncation error in 35 is of order Axnz as we already found for the explicit and fully implicit cases with the exception of a special case when 9 05 the midpoint method for which the local truncation error is of order Ax For an ODE it is possible to write a general code representing 35 with 9 simply an input parameter In this way one can leave it to the user to choose which Euler method to use Performance of the method varies with the equation function f being solved the step size Axquot and the choice of weight 9 For example we ll see later that for 05 5 9 5 1 the Euler method is New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology unconditionally stable but that for 0 5 9 lt 05 stability is conditional constrains Axquot lt some critical value in order not to blow up Consider the linear rstorder ODE p 26 of text y x px y rx where p and r are functions depending on the independent variable x Then f r 7 py and equation 35 becomes y s y Ax1 6w py6 r py 36 Collecting terms involving yn1 which is implicit in 36 I1Axnt9p1yn1E1Axn1t9pnynAxn1t9rnt9rnll 37 The rst term on the RHS represents the response to the IC in an initial value problem while the second term represents the effects of the forcing r mimicking the exact solution see 5 on p 78 of these notes or p 28 of the text Continuing to simplify for yn we get the algebraic equation l Ax 1 6 l t9rt9r ym E myn WC HM 9p HM 9pm 38 What does this have to do with solving partial differential equations PDEs This general Euler approach is commonly used to solve the time derivative in PDE solutions of groundwater ow heat transport and solute mass transport In those cases the space dimensions are usually handled by nite difference or nite element methods to be introduced later In the Euler approach to ODEs we divide by the coef cient 1 Axquot t9 pm1 on the lefthand side of 37 in order to get 38 In the Euler approach to the time derivative of PDEs the lefthandside coefficient is actually an mxm matrix where m is the number of spatial node points If 9 0 this coefficient is the identity matrix I Since 1391 I this division operation is trivial for the forward method For 9 gt 0 the coef cient is more complex and simultaneous equation solvers like the Gauss Elimination are needed We ll introduce these solvers a little later when we start to deal with 2quotd order ODEs and boundary value problems BVPs The choice of 9 then effects how you solve the spatial part of the PDE Many PDE codes allow fully implicit and midpoint solutions for the time derivative Some allow any value of 9 2 05 or even 9 gt 0 In these codes programming both the option for an explicit solution 9 0 and the option for implicit solutions 9 gt 0 is cumbersome since only the implicit solutions require simultaneous equation solvers for space Consequently that combination is unusual In some hydrology codes different solvers are used in different parts of the spatial domain So called mixed implicitexplicit time derivatives are used These codes determine where in space the explicit solution 9 0 is stable and use that scheme in that part of the domain an implicit solution say 9 l is used elsewhere This allocation is calculated dynamically Midpoint Method In the midpoint method 9 05 and 34 and 35 become 103 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology ym y 1 E 39 Ax 2 f f and NH y 2y 2 f f 40 This method is said to have less local truncation error than Euler methods with other choices of 0 and the same as the improved or modi ed Euler method Huen s method Why The exact ODE at the midpoint is y xn12fxn1Z 41 Approximate y xn1z by 30 with 9 05 l rm 2 2Lv39x y39x 42 and replace the derivatives by their series expansions 6 and 15 keeping all the terms Axl y xn quot N 1 yxn1yxn Ax y xn12 2 Ax 2 3 y xn H l l MH2 2 Mn 2 y 9 3 Notice that the rst terms of each expression on the RHS is the same The second terms are very similar but with a change in sign They would cancel out if the second derivatives were the same at the two endpoints Assume that they are approximately equal and neglect these terms The third terms have the same sign and will be additive And it continues to alternate like this with every other term approximately zero The result is 1yx yxg ywxmi 43 yxn1yxnAX nyvllg 44 y xn12 E 3 n When we truncate this exact expression and substitute it and 32 into the original equation 41 the numerical difference and the solution for yn1 are higher order yn1 yn l 2 TX 2 f f 0w 45 y s y f f 0Ax3 47 That is the midpointmethod local truncation order for yn1 is Aan while it is of order Axnz for all other Euler methods For this reason the midpoint method is often programmed into code and 104 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology is preferred by many users While theoretically it might have a lower truncation error actual performance is usually not so obviously better Stability of Euler Methods Let s examine stability by focusing on the linear ODE p 27 of text y x px y rx where p and r are independent variables It will be easier to focus if we set the forcing r to zero and if the coefficient p is a constant That is consider the linear ODE y x p y 0 48 This is a version of 1 withf py The exact solution of 48 is y C6 an exponential response If we take the initial condition IC as yxn and predict yxn1 we can evaluate the constant C and write the particular solution yxn1 yxn 6Wquot 49 where Axquot xquot 1 7 xquot The dependent variable decays exponentially in value from the IC at xquot to the nal condition at xn The general Euler method approximate solution for the same situation is given in 38 with r0 and pconstant or ymyn ngmn yA 50 where A is a socalled ampli cation factor that is attempting to mimic or approximate the exponential behavior in 49 1 1 9PAxn 1 19 P Axquot A 51 How well doesA mimic e39m First note that 0 Se PAX S 1 while these bounds don t apply to A Then compare plots of A and e PAX for various values on say 9 0 025 05 075 10 using Matlab For which values on and pr is lAlgtl For which 9 does the error grow faster than the solution decays stiffnesss and becomes unstable small perturbations grow Under what conditions is the solution always unconditionally stable Here is a sketch 6 See p 896 of text 105 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology exponential response 1910 Kl as A gtoo j stability limit 00 If lAlgtl we have uncontrolled growth and instability Why Recall that yn1 z ynA thus yl z yoA where yo is the IC yZ ylA ZyoAz y3 Z y0A3 etc so thatynZ yoAquot If lAlgtl then Aquot will blow up and with it y We want to mimic decay not growth We need stability That requires that lAlltl In the diagram note the values of A as pr approaches in nity Call these A00 From the diagram if 0 5 9 lt 05 stability depends on the value of pr The method is conditionally stable leoPl If 05 5 i9 5 l the method is unconditionally stable le451 The midpoint method is borderline conditionally stable leoll In practice it is desirable to choose 05 5 i9 5 1 so solution never blows up 9 05 and l are the most popular choices Restricted step size for forward method What if you choose to use the explicit or forward method with 9 0 From the sketch you can see that if prn gt 2 then the method will blow up This restricts the step size Axquot so that prn 5 2 0r Axquot 5 2p to avoid instabilities Recall from our discussion of the analytical solution for this problem that the efold drop decay constant or time constant is Up ie e39p Axquot e39A quot d ay quotsquot quot39 Thus Axquot 5 2 X decay constant We want the step size to be smaller than the decay constant anyway so this doesn t seem so onerous despite what our textbook says However when solving the time domain for big hydrology PDEs see box on p 103 of these notesthe time constan is not that of the system but of the spatial nite difference grid block or nite element used to handle the spatial part of the problem In that case our p is equivalent to the eigenvalue of that block or element in the conductivity coef cient matrix and strongly restricts the solution We ll visit this later 106 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology RungeKutta Methods Background text pp 892898 85 Gilat and Subramaniam 2008 Once again we are solving the IV problem y fxy yxoy0 1 without restricting it to linear equations Recall the general approach for singlestep explicit methods xnl m xn Axquot l6ab ym e y AX slope The explicit Euler method used the slope of y at xquot The RungeKutta RK method calculates the slope from a weighted average of estimates of the slope of y at several points within the interval xquot to xquot 1 Higher order RK methods use more points and have smaller truncation error For example the socalled classical RK method is of 43911 order globally and uses four points This is the method described in our textbook The second order RK method uses only two points and is 2quotd order globally Recall that the backward and forward Euler methods are 1st order globally The modi ed Euler method and the midpoint Euler method are 2quotd order globally just like the 2quotd order RK method The difference between the methods is in the location of the points within the interval that are used to calculate the slope Second order RK methods not covered in the textbook The general form is given by yn1 yn crk1 czkz Axn with k1 xnyn and k2 xnazAxn ynb21k1Axn and where Cl 02 a2 and bu are constants The modi ed Euler method is a special case of the 2quotd order RK scheme particular values of these constants with both cl 12 02 12 a2 l and bu 0 Third order RK methods not covered in the textbook The general form is given by yn1 yn crk1 czkz csks Axn with particular de nitions for the k s and eight constants to assign values to The classical 3rd order RK method has particular values of the constants and is related to Simpson s rule of numerical integration 85 Gilat and Subramaniam 2008 107 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Fourthorder RK methods pp 892894 of textbook The general form is given by yn1 yn Clk1 Czkz Csks C4k4 Axn with k1 xnyyn k2 xna2Axny ynb21k1Axn k3 xna3 AXm ynb31k1AXn b3zk2AXn k4 xna4Axny yn b41k1Axn b42szXn b43k3Axn In the classical 4th order RK method 01 CZ C3 26 C4 02 12 a3 12 a4 1 and 521 12 b31 0 b32 12 541 0 542 0 b431 From the text an explanation and pseudocode A method of great practical importance and much greater accuracy than that of the improved Euler method is the Clusxit39al Rllllgt iKllITU NWT101 of mrrll order which we call briefly the Runge Kutta method1 it is shown in Table 214 We see that in each step we first compute four auxiliary quantities k1 kg 3 4 and then the new value 1 The method is well suited to the computer because it needs no special starting procedure makes light demand on storage and repeatedly uses the same straightforward computational procedure It is numerically stablei Note that if f depends only on X this method reduces to Simpson s rule of integration Sec 195 Note further that M k4 depend on I and generally change from step to step Table 214 Classical Runge Kutta Method of Fourth Order ALGORITHM RUNGEiKUTTA f x0 0 II N This algorithm computes the solution of the initial value problem 7 ftx y 39 0 yo at equidistant points x l x o 11x2 x o 21 x N x0 Nb here f is such that this problem has a unique solution on the interval 360 rN see Sec 17 INPUT Function f initial values X0 yo step size I1 number of steps N OUTPUT Approximation 39H1 to the solution 39rn1 at rn A39O n lI where1 O l 39 N 108 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Fern 01N1do k1 hfrn N k2 hftrn 1 kl k3 I1frn ll y kz k4 hfrn II yn k3 rnH rn I1 57 yn k1 2k2 2k3 k4 OUTPUT rn n1 End Stop End RUNGEKUTTA Classical Runge Kutta Method Apply the RangeiKutta method to the initial Value problem 4 in Example 1 choosing I1 02 as before and computing five steps Solution For the present problem we have x y 39 r Hence k1 02xn jun k2 021 01 1 05k1 k3 02xn 01 yn 05k2 k4 021 02 1 k3 Table 215 shows the results and their errors which are smaller by factors 103 and 10 1 than those for the two Euler methods See also Table 216 We mention in passing that since the present k1 r 394 are simple operations were saved by substituting k1 into 2 then k2 into k3 etc the resulting formula is shown in Column 4 of Table 215 Table 215 RungeKutta Method Applied to 4 n X V 022140L 31 Exact Values 6D 106 X Error 39 1 00214 y e 7 r 7 l of 0 00 0 0021 400 0000 000 0 1 02 0021 400 0070 418 0021 403 3 2 0 4 0091818 0130 289 0091825 7 3 06 0222 107 0203 414 0222 1 19 12 4 08 0425 521 0292 730 0425 541 20 5 10 0718 251 0718 282 31 Comparison from the text YIN 116 Comparison of the Accuracy of the Three Methods Under Consideration in the Case of the Initial Value Problem 4 with h 01 Error tquot quot quotI quot l Euler improved Euler Rungc Kulla Table 211 Tahh 213 Table 215 02 1021 403 0021 00014 0000 003 04 0091 2415 0052 00034 0000 007 06 0223 11 00 11 00015 0000 011 03 0415 541 0152 00101 0000020 10 0713 2113 0229 010156 11000031 109 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology RungeKuttaFelhberg RKF Error and Step Size control with RK The idea of adaptive integration Sec 195 has analogs for RungeiKutta and other methods In Table 214 for RK RungeiKutta if we compute in each step approximations and with step sizes I1 and 211 respectively the latter has error per step equal to 25 Z 32 times that of the former however since we have only half as many steps for 211 the actual factor is 252 16 so that say 52111 x 16601 and thus Xli 7 217 E217 7 011 g 16 7 l h1 Hence the error 5 EU for step size It is about 10 e 7 i where e 1 7 ya as said before Table 217 illustrates 10 for the initial value problem 11 y y 71 7 n2 2 10 1 the step size I1 01 and 0 E 39 E 04 We see that the estimate is close to the actual error This method of error estimation is simple but may be unstable Table 217 Runge Kutta Method Applied to the Initial Value Problem 11 and Error Estimate 10 Exact Solution y tan x x 39I r39 Error Actual Exact Step size II Step size 2h Estimate 10 Error Solution 9D lt2 00 1000 000 100 1000 000 000 0000 000 000 0000 000 000 1000 000 000 01 1200 33 1 589 0000 000 083 1200 334 672 02 1402 709 878 1402 707 408 0000000 I65 0000000 157 1402 710 036 03 1609 336 039 0000000 210 1609 336 250 04 1822 792 993 1822 788 993 0000 000 267 0000000 226 1822 793 219 One can run RK with these two different step sizes compare the results and if the error estimate is too large reduce step size and do it again In application this is done adaptively at each step The RungeKuttaFelhberg RKF algorithm uses a more sophisticated and ef cient approach to stepsize control It uses two RK methods of different orders for each step and picks those methods to minimize the number of arithmetic operations If the error estimate is too large the step size is adjusted The method has a global truncation error of 5 11 order The text covers this pp 894895 because RKF is very popular eg in Maple and Matlab In hydrology it is the second most popular way to track particles in path line and travel time simulations just after Euler methods llO New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Multistep Methods see text 213 Single step methods use only one previous step at xquot in order to predict the value of the dependent variable yn1 at the next step xn1 Multistep methods look further backwards to solutions and information at older steps xn1 xmz etc Recalling that we have certain 39 on the J r J variable y this provides more information about how y is behaving and a more accurate prediction of yn1 Different multistep methods use different amounts of older information number of previous steps used and employ that information differently Because they use more past information each prediction step is more computationally expensive more arithmetic operations and storage The text reviews the AdamsBashforth and AdamsMoulton multistep methods both shown for their respective 4Lh order versions One of problems with implementing these methods is getting them started At the beginning of a computation the rst initial value IV yo is at x0 with right side f0There are old values of y and f at x1 xz etc Consequently you have to start the solution with another method say RK or Euler until you are far enough into the domain to have all the old values that you need with the multistep method Multistep methods are becoming more common in applied codes so you should be familiar with the basic concepts but there is no need to become expert in this introductory course New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Part2 Ordinary Differential Equations ODEs This is new material see Kreyszig Chapters 16 and related numerics in Chaps 19 20 201207 and 211213 Fundamentals of Differential Equations The calculus problems we ve reviewed have mostly been involved with nding the numerical value of one magnitude or another For example we ve sought the location x0 and value xo of the maxima of a function the locations x where x0 roots of fx the slope derivative of fx at x or the value of a de nite integral 17 f xdx But now we are going to solve engineering and science problems where the unknown is itself a function that is we are going to develop and solve models that express the dependence of certain variables on others Petrovskii 1956 What kind of questions can these models answer Here are some examples from hydrology How will drawdown in an aquifer change with distance from a pumping well How will the discharge from a lake owing over a weir vary in time How does the velocity profile in the atmospheric boundary layer vary with distance above the ground surface How does moisture content vary with depth In these cases the dependent unknown variables are respectfully drawdown L1 discharge L3 velocity Lt and moisture content L3 L3 The respective independent known variables are distance L time t elevation L and depth L The problem of nding these functions is most often addressed by solving difkrential equations that is equations in which not only the unknown functions occurs but also its derivatives of various orders Examples L ahahoP alt d2 1 Qv 2 W dxz Ni 1 T U l C l I In the rst example the unknown is denoted by the letter h which represents the water level in a lake and the independent variable by the letter t which represents time As there is one independent variable it is an ordinary differential equation ODE As there is only one rst order derivative it 1The notation L for length and its like eg M for mass or t for time is used to indicate units or dimensions of a variable similar to using standard units like In meter 5 second or kg kilogram So velocity units can be indicated by Lt or ms39l Among other things we use units to keep track of dimensional homogeneity such that units balance in all terms on each side of an equation See next page New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology is a first order equation see below The other symbols a ho P2 denote parameters with no dependence on the unknown h Mathematically its possible that all three vary with the independent variable t but in hydrologic practice only no and P typically vary with time In the second example the unknown is denoted by the Greek letter representing hydraulic head in a con ned aquifer and the independent variable by the letter x representing location As there is still just one independent variable it is also an ordinary di erential equation ODE However as a second derivative of is the highest order derivative it is a second order equation The symbols A 0 and Qw denote parameters that are known In the third example the unknown is denoted by the letter C representing solute concentration in a owing river and there are two independent variables denoted by the letters x space and t time As there are now two or more independent variables it is a partial differential equation PDE The highest order derivative in x de nes it as a second order equation in space while the highest order derivative in tde nes it as a first order equation in time Later we ll see that this particular type of PDE is called mixed hyperbolicparabolic PDE The symbols v and D denote parameters the mean river velocity and a dispersion coef cien In each equation the left side contains the unknown and its derivatives while the right side contains known loadings or forcings For example in the rst equation the forcing is precipitation P Lt In the second equation it is pumping Qw Lt and head in a vertically adjacent aquifer 0 L There is no forcing in the third equation later we ll see that forcing for that equation comes from the boundary and initial conditions instead We often rewrite different equations so that they t this template unknowns and their derivatives on the left forcings on the right As we shall see below mathematicians take this convention a step further In short the first equation is a first order ODE the second equation is a second order ODE and the third equation is a second order in space 7 first order in time PDE In this section we concentration only on problems with one independent variable that is on ODEs Dimensional homogeneity Each of these DEs is and must be dimensionally homogeneous We can take advantage of this to help detect errors and to improve solutions For example the second equation is 2 V 1 i i o QW x l l T Here head and distance x have units of length L parameter l L2 includes TB K the in uence of aquifer transmissivity T Lzt and aquitard leakance B K LtL where B is aquitard thickness L andK is aquitard vertical conductivityLt and Qw Lt is pumping per unit area All terms in the second equation have units of lL Order of equation text 1l An ordinary differential equation of order n in one unknown function y is a relation of the form leyxy39 xy xyquot 96 0 Between the one independent variable x and the quantitites 2 The physical meaning of parameters like these will become known to you later 59 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology yx y39 x y39 39 XLm yquot x The order of the ODE is the order to the highest derivatve of the unknown function appearing in the differential equation We ll focus on 1st and 2quotd order ODEs respectively Chaps 1 and 2 of the textbook In our applications time is usually the independent variable in lSt order ODEs while it is almost always space in 2quotd order ODEs The major exception to these observations is in geophysics involving seismic electromagnetic and other signals where 2quotd order ODEs in time are also encountered Operator We can write an equation in the form of an operator L where Lunknown forcing or L forcing In the examples above the operators are d 1 L7a d2 2L dx 6 6 62 3 L v DES An operator can be an algebraic differential or integral operator Our examples only demonstrate differential operators Linear dependence There are tremendous advantages of a linear model in deriving and applying solutions For example if we have a linear model then we can scale or convolute the model result to represent a new forcing that is a new solution without having to resolve the equation If we have a linear model we can use certain solution methods that take advantage of linearity such as LaPlace or Fourier transform methods Even if the model is nonlinear if the nonlinearity is mild we might be able to linearize approximate the model and still take advantage of linearity On the other hand there are certain problems that have a suf cient nonlinearity that it must be addressed directly Consider chaotic behavior as in climate models which is the result of one kind of nonlinearity that cannot be ignored In vadose zone hydrology the multiphase ow equations are highly nonlinear The landsurface energy balance involves highly nonlinear radiation terms Taking advantage of the operator notation there are three conditions necessary for an operator and its differential equation to be linear Letting x and y represent unknowns and a and b constants these are a Lax aLx New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology b Lx y L06 My 0 Lax by aLx bLy As you can see the last condition covers the other two Let s take a look at two algebraic examples Example a bxc 0 has the operator L b so that Lx c Show that L meets the requirements for a linear operator Example b ax2 bxc 0 has the operator L a 2 b so that Lx c Show thatL does not meets the requirements for a linear operator That is this is a nonlinear algebraic equation Aside if b0 show that the operator remains nonlinear in x but that a transform of variables to zx2 leads to a new problem azc 0 which is linear in z We often look for such simple transformations to exactly linearize a problem such transformations are common in vadose zone hydrology and in ltration models used for hillslope or watershed hydrology and for simple models of phreatic aquifers through the socalled linearized Dupuit Approximation Now let s examine operators for differential equations Example 0 ah aho P has the operator L a If dis a constant is this operator linear Example 0 ah P has the operator L a 1 If a is a constant is this operator linear 2 Example e What conditions on v and D are needed for L g v D x linear operator Hint regarding linearity look to see if the unknown is of the same form order of magnitude power in each term but formally you need to apply rule 0 above Concept of a solution text 11 Given the rst order ode F x yx y39 x 0 4 a function y hx is called a solution to 4 on some open interval a lt x lt b including ab oo 00 if 61 to bea New Mexico Tech Hydrology Program Hyd 510 Quantitative Methods in Hydrology hx is de ned and differentiable throughout the interval and if the equation becomes an identity if y and y are replaced with h and h respectively The curve graph of hx is called a solution curve Verify a solution Always test a solution by substituting it into the original ode and solving Solution curves text l l l The ODE Iquot dyili c0s can be solved directly by Integration on both sides Indeed using calculus we obtain i39 f ms r il39 sin 39 t39 where t39 is an arbitrary constant This is aftwily of solutions Each Value of t for instance 275 or 0 or 78 gives one of these curves Figure 2 shows some of them for t39 i l 0 2 3 4 The family of solutions is called the general solution see next page Each of these many solutions is equally valid How do we choose a single one of these solutions from within family Such a single solution is called the particular solution We need some additional information In this example of a 1st order equation we need just need one value of the solution yx0 at some location within the domain ice to nd the particular solution In fact we almost always pic alocation on the ends ofthe domain For 1st order ODEs this is the initial location x0 and the value yx0 is the initial value of the solution yo yr1392 39 Fig Solutions y sin x of the ODE y39 05 x In the example to the right the selected solution is the one shown in the highlighted line with a particular value yx0 y7I 2 used to select the particular highlighted solution 1quot Order ODEs Textbook Chap 1 Exponential rowth exponential decay text 11 From calculus we know that 39 03 t39 any constant has the derivative chain rule 39 341quot 3r e l q r This shows that V39 is a solution ol r J L Hence this ODE can model exponential growth tor instance 01 animal populations or coloni s ot bacteria It also applies to humans for small populations in a large country leg the United States in early tiniest and is then known isillriltliiix s Itiwl We shall say more about this topic in Sec 1 i Siiii39il I o V 01V With a minus on the right has the solution 7 u 2 Hence this ODE niotlels exponential decay for instance of a ratlinacliye substance see Example 5 Figure 3 sliOVx s solutions for some positch Can you find what the solutions look like for negative 1quot 62 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology O 2 4 6 8 10 12 14 t Fig 3 Solutions ofy 702v in Example 3 General solution arbitrary constant c and the family of solution curves as in Figures 2 amp 3 Particular solution a particular value ofc and a particular curve eg one of the curves in Figures 2 amp Introduction to initial value problems text 11 Generally a particular solution to a first order problem is obtained by specifying an initial condition 1C yxo yo a known value of y at xxo We then seek the solution y for xgtxo That is we use the IC to determine the particular solution This type of problem is called an initial value problem y39 fxy yxo y0 5 which includes the ODE left and IC right Since we are solving a first order equation we need only one IC Later we ll see that a second order ODE requires two ICs for an initial value problem Geometric Meaning of y fxy text 12 A first order ODE y39 f x y 1 has a simple geometric interpretation A solution curve of 1 that pass through a point x0y0 must have at that point the slope y39 x0 equal to the value offat that point that is y39xo fx0gtyo You can indicate the direction of solution curves of l by drawing short straightline segments lineal elements in the xyplane see figure and then fitting approximate solution curves through the resulting direction field New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology r rrr oclines of xy VKgik const alon hich the slope is fiant a By Matlab or by a CA8 b By isocllmes computer algebra system like able r Fig 7 Direction field ofy xy Why do this especially since it is of limited accuracy 1 You don t have to solve l which is especially convenient when the solution to l is challenging 2 You see in graphical form the entire family of solutions and their propertieses New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Separable ODES previously mentioned in our calculus review on pp 3537 of these notes Many practically useful ODEs can be reduced to the form D wW m by purely algebraic manipulations Then we can integrate on both sides with respect to X obtaining 2 gm quot Ix fftx Ix c On the left we can switch to y as the variable of integration By calculus quot le cl so that m mm mma This method of solution is called the method of separation of variables and l is called a separable equation x appears only on the right and y appears only on the left Example The ODE y l 2 is separable because it can be written as W lyZ By integration the solution is arctany x e or y tanx e New Mexico Tech Hyd 51 Hydrology Program Quantitative Methods in Hydrology Modeling Mixing Problem Mixing problems occur quite frequently in chemical industry We explain here how to sole the basic model inol39ing a single tank The tank in Fig 9 contains 000 gal of water in which initially 100 lb ofsalt is dissohed Brine runs in at a rate of IO galmin and each gallon contains 5 lb of lissoyed salt The mixture in the tank is kept uniform by stirring Brine runs out at IO galmin Find the amount of salt in the tank at any time I Solution Step 1 Setting up a model Let rm denote the amount of salt in the tank at time I Its time rate of change is Iv Salt inflow rate 7 Salt outflow rate Balance lawquot 5 lb times I0 gal gives an inflow of 50 lb of salt Now the outflow is It gal of brine This is I0I000 00I 1 ofthe total brine content in the tank hence 001 of the salt content t that is 00lrr Tlnls the model is the ODE 4 y 7 50 7 0013 7 700107 5000 Step 2 Solution oftw model The ODE 4 is separable Separation integration and taking exponents on both sides gives dv 739 7 7001 m In Ir 7 5000i 7 7001 11 r 7 5000 7 two01 r e 3000 Initially the tank contains I00 lb of salt Hence j390 I00 is the initial condition that will give the unique solution Substitutingr 100 and I 0 in the last equation gives I00 7 5000 L39 Hence L 74900 Hence the amount of salt in the tank at time I is 5 vol 7 5000 7 4900n 001 This function shows an exponential approach to the limit 5000 lb see Fig 9 Can you explain physically that rl should increase with time 7 That its limit is 5000 lb 7 Can you see the limit directly from the OD 7 The model discussed becomes more realistic in problems on pollutants in lakes see Problem Set 15 Prob 27 or drugs in organs These types of problems are more dif cult because the mixing may be imperfect and the flow rates in and out may be different and known only very roughly 3 5000 4000 3000 2000 1000 loog l l l l 0 100 200 300 400 500 t Salt content ytt Tank Fig 9 Mixing problem in Example 3 55 New Mexico Tech Hyd 51 Hydrology Program Quantitative Methods in Hydrology Radiocarbon DatingZ In September 1991 the famous Iceman Oetzi a mummy from the Neolithic period of the Stone Age found in the ice of the Oetztal Alps hence the name Oetzi in Southem Tyrolia near the Arrstrianrltalian border caused a scientific sensation When did Oetzi approximately live and die if the ratio of carbon 6C to carbon SC12 in this mummy is 525 of that of a living organism Physical Information In the atmosphere and in living organisms the ratio 0139 radioactive carbon 6C made radioactive by cosmic rays to ordinary carbon 6c is constant When an organism dies its absorption of C by breathing and eating terminates quot quot 39 we nFa Fossil by comparing 39 ratio in the fossil with that in the atmosphere To do this one needs to know the halflife of 5C which is 5715 years CRC Handbook of Chemistry and Physics 83rd ed Boca Raton CRC Press 2002 page 11752 line 9 Solution Modeling Radioactive decay is governed by the ODE y ky see Sec 11 Example 5 By separation and integration where t is time and yo is the initial ratio of SC to SC l 7y k 111 In Iyl U c y yoek39r yo 3 Next we use the halflife H 5715 to determine k When 1 H half of the original substance is still present Thus ln 05 0693 kl an yoe 05m e 05 k H 5715 0000l213 Finally we use the ratio 525 for determining the time t when Oetzi died actually was killed 111 0525 let 70 0001213 t 1 A A e e 0 525 t 410001213 5312 mwer About 5300 years ago Keep in mind for this and all All of these examples have hydrologic relevance As you can see ho mlogou v I orderPrOblemquot from them many hydrologic problems have an exponential response to some IC either decay or growth Rate Constant 0693 X half39hfe It 0693 x H Heating an Office Building Newton s Law of Cooling Suppose that in Winter the daytime temperature in a certain office building is maintained at 70 F The heating is shut off at 10 FM and turned on again at 6 AM On a certain day the temperature inside the building at 2 AM was found to be 65 F The outside temperature was 50 F at 10 PM and had dropped to 40quotF by 6 AM What was the temperature inside the building when the heat was turned on at 6 AiM 7 Physical information Experiments show that the time rate of change of the temperature T of a body B which conducts heat well as for example a copper ball does is proportional to the difference between T and the temperature of the surrounding medium Newton s law of con ing Solution Step 1 Setting up 11 model Let Tt he the temperature inside the building and TA the outside temperature assumed to be constant in Newton s law Then by Newton s law dT 6 I MT 7 TA Such experimental laws are derived under idealized assumptions that rarely hold exactly However even if a model seems to t the reality only poorly a in 39 itma still iv 39 39 39 39 To see how good a model is the engineer Will collect experimental data and compare them with calculations from the model New Mexico Tech Hyd 510 Step 2 General solution We cannot solve 6 because we do not know TA just that it varied between SO F and 40 F so we follow the Golden Rule fyau cannot solve your problem try to solve a simpler one We solve 6 with the unknown function TA replaced with the average of the two known values or 45 F For physical reasons we may expect that this will give us a reasonable approximate value of Tin the building at 6 AM For constant TA 45 or any other constant value the ODE 6 is separable Separation integration and taking exponents gives the general solution dT T 45 k It ln T 7 45 k1 at T0 45 ce 6 e Step 3 Particular solution We choose 10 RM to be I 0 Then the given initial condition is T0 70 and yields a particular solution call it I p By substitution 70 45 cequot 70 c 70 45 25 TF0 45 258 Step 4 Determination ofk We use T4 65 where l 4 is 2 AM Solving algebraically for k and inserting k into TF0 gives Fig 10 1104 45 2551 65 e4 t 08 k 11 In 08 70056 rpm 45 25er0 563 Step 5 Answer and interpretation 6 AMA is r 8 namely 8 hours after 0 PM and Most of these problems are forward problems Given the ODE 1C and parameter values how does the unknown vary with the T 8 45 25 56 39 8 61 F 1393 Hence the temperature in the building dropped 9 F a result that looks reasonable 5 independent variable 7 7o 7 Step 4 is an example of an 68 r invase orp arameter estimation problem Given the ODE 1C an observatioions of the dependent variable what is an estimate for one or more ameters E7 2 4 6 8t Fig10i Particular solution temperature in Example 4 Leaking Tank Outflow of Water Through a Hole Torricelli s Law This is another prototype engineering problem that leads to an ODE It concerns the out ow of water from a cylindrical tank with a ho e at the bottom Fig 11 You are asked to find the height of the water in the tank at any time if the tank has diameter 2 m the hole has diameter 1 cm and the initial height of the water when the hole is opened is 225 in When will the tank be empty Physical information Under the influence of gravity the outflowing water has velocity 7 ut 0600V2ght Toi ricelli s law where ht is the height of the water above the hole at time t and g 980 cmsecz 3217 ftsee2 is the acceleration of gravity at the surface 0139 the eart Solution Step 1 Setting up the model To get an equation we relate the decrease in water level 121 to the out ow The volume AV of the outflow during a short time At is AV Av At A Area of hole Borda s coe icient ofcontraction C06U so thatAC7rDz4 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology AV must equal the change AW 0139 the volume of the water in the tank Now AV 3 Ah B Crosssectional area of tank where Az gt 0 is the decrease of the height 11 of the water The minus sign appears because the volume of i the water in the tank decreases Equating AV and AW gives 3 Ah Av At We now express 1 according to Torricelli s aw and then let At the length of the time interval considerad approach O this is a standard way of obtaining an ODE as a model That is we have An A A 7 A 7 x Ar B I B 0600 2ght and by letting A gt 0 we obtain the ODE l A What kind of ode is this 7 17 2656 E W Order linearity othe Where 2656 0600 2 980 This is our model a rstorder ODE Step 2 General solution Our ODE is separable All is constant Separation and integration gives dh A 7 2656 7 dt and 2f 2 W B Dividing by 2 and squaring gives it c M 1328AtB2 Inserting 1328MB 1328 05211100271 0000332 yields the general solution 2656 A B t 112 6 00003322 Step 3 Particular solution The initial height the initial condition is 10 225 cm Substitution of t 0 and h 225 gives from the general solution 02 2 c 1500 and thus the particular solution Fig 11 hpo 1500 7 000033202 Step 4 Tank empty hpt 0 il39t HBO0000332 45 181 sec 126 hours Here you see distinctly the importance of the choice of unitsiwc have been working with the Cgs system in which time is measured in secendsl We used g 980 crnsecz Step 5 Chec Check the result I How does the velocity v and discharge VA vary with time 0 10000 3000 50000 2 Tank Water level Ml in tank Fig 11 Example 5 Outflow from a cylindrical tank leaking tankquot Torricelli s law New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Reduction to Separable Form Certain nonseparable ODEs can be made separable by transformations that introduce for y a new unknown function We discuss this technique for a class of ODEs of practical importance namely for equations 1 V 3 3 f X Here f is any differentiable function of r such as sin r rX4 and so 011 Such an ODE is sometimes called a unnogeneom ODE a term we shall not use but reserve for a more important purpose in Sec 15 The form of such an ODE suggests that we set 7 It thus 9 m and b Jroduct differentiation r u x L1 YI Substitution into quot f t then gives u u u or u x u 7 it We see that this can be separated It I 10 7 f Lt u r Exact ODEs Integrating Factors Basics of Exact ODEs A firstorder ODE M y N quot 0 written as use ly Itquot l as in Sec 1 2 1 M y Ix Nt y Iy O is called an exact differential equation if the differential form M39 y lt39 N39 y zly is exact that is this form is the differential Em 39 t 2 IL Ix i Ir rl Y 39 of some function 11 y Then I can be written In 0 New MexicoTech Hyd 510 Hydrology Program Quantitative Methods in Hydrology By integration we immediately obtain the general solution of l in the form 3 u y v The function u needs to satisfy 61 7 611 4ab a 7M a N The necessary and suf cient condition for exactness is see text p 20 8M 0N EI ax 39 5 This condition is not only necessary but also sufficient for 1 to be an exact differential equation We shall prove this in Sec 102 in another context Some calculus books eg Ref GRl 1 also contain a proof If 1 is exact the function 14 y can be found by inspection or in the following systematic way From 4a we have by integration with respect to x 6 LI chx k in this integration is to be regarded as a constant and k plays the role of a constant of integration To determine k we derive Flut from 6 use 4b to get quotf17 and integrate tlkd to get k We can equally apply 4b to get the following If 6 doesn t work then try 7 7 u j Ndy ltx See text page 23 for example New Mexico Tech Hydrology Program Quantitative Methods in Hydrology Reduction to Exact Form Integrating Factors We multiply a given nonexact equation 12 Ptx y be QM 7 17 O by a function F that in general will be a function of both x and y We want the result to be an exact equation 13 FPdX FQ d 0 so we can solve it as just discussed Such a function F x y is then called an integrating factor of 12 See text 1 24 for example How to Find Integrating Factors 1 l6 Z R where R E dF 15 8Q iv 39Lr Integrating Factor Fx 1f12 is A39LICI 1 qu TIz rigII ride R 37 l6 depends only on x39 then 12 has ill integrating fm39mr F F39 whim ix obtained by illk gl39ulilig 16 and taking exponents 7n mti sides 17 Fr exp J RU d39 Integrating Factor Fy If 12 is YLICl that the rign side Ri 0f18 depends 0in on y Then 12 lay an integmring icmr F F y which is 0Rlillt t 0lil 18 in Iiiefm m 19 F exp R dy New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 2 Matrices Matrices matrix a rectangular array of number A matrix is usually identified by a boldface capital letter The size of a matrix is m x n read m by n where m is the number of rows in the matrix and n is the number of columns An entry in a matrix is denoted by a lowercase letter with two subscripts indicating the row and column of the entry eg an is the entry in the ith row and jth column of matrix A Matrices are used to solve systems of linear equations General matrix 11 12 17 A 21 22 I 2 aml amZ am Example A is a 2x3 matrix A 1 5 32 5 r 1111 112 1332 1217 122 2 23TC 7 J2 7 Special Matrices Square matrix a matrix that has the same number of rows and columns ie n x n 4 1 0 2 1 A B 1 2 2 4 5 3 4 5 Diagonal matrix a square matrix whose only nonzero entries are on the diagonal ie an 2 O fori j 0 0 0 A OONO p 1000 3 0 0 0 Tridiagonal matrix a square matrix whose only nonzero entries are on the diagonal andjust above andjust below the diagonal Hyd 510 Quantitative Methods in Hydrology Lab Notes 2 Matrices New Mexico Tech Hydrology Program 3 2 0 0 0 1 5 1 0 0 A 0 4 4 3 0 0 0 1 7 9 0 0 0 2 5 Upper Triangular Matrix U a square matrix that has all zeroes below the diagonal Lower Triangular Matrix L a square matrix that has all zeroes above the diagonal 2 1 1 0 2 0 0 0 0 3 7 2 1 5 0 0 U L 0 0 6 1 3 7 2 0 0 0 0 5 4 3 0 8 Identity matrix I a diagonal matrix with all ones on the diagonal 1 0 0 0 I0100 0010 0001 Symmetric matrix a square matrix with an 11 1 2 1 A 2 3 0 1 0 7 Matrix Operations Let A 11 12 13 B b11 712 I713 21 22 23 721 722 723 C C C 11 12 d11 d12 C21 022 D 21 22 C31 C32 Matrix addition two matrices of the same size dimensions can be added by adding the entries in the same positions If Z AB then zij aijbij 11 1711 12 712 13 713 Z A B 21 1721 22 1722 23 1723 Matrix subtraction two matrices of the same size can be subtracted by subtracting the entries in the same positions If Z AB then Z Iiibu New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 2 Matrices 11 711 12 712 13 713 b 21 21 22 b22 23 b23 ZA B Scalar multiplication each entry of the matrix is multiplied by a scalar If Z 01A then Zijaaij Z aA W11 W12 W13 W21 W22 W23 7 Matrix multiplication an m x n matrix can be multiplied by an n xp matrix If ZAC then zi 2a k ck for k1 i12mandj12p Z AC 11c11 12021 13031 11012 12022 13032 21011 22021 23031 21012 22022 23032 011 11 c12 21 011 12 c12 22 011 13 c12 23 Z CA 021 11 c22 21 021 12 c22 22 021 13 0226123 031 11 c32 21 031 12 c32 22 031 13 0326123 Note that the order of multiplication is important ie AC is not necessarily equal to CA The number of columns of the first matrix must be equal to the number of rows of the second matrix in order to multiply the two matrices If A isnxmandBismxpthenZABisnxp Matrix transpose If AT is the transpose of A then the rows of A are the columns of AT 11 21 T A 12 22 a a 13 23 Note that if A is symmetric then ATA Also note that the transpose of an upper triangular matrix is a lower triangular matrix and the transpose of a diagonal matrix is diagonal Determinant the determinant of a square 2x2 matrix is lDl 11d22 12d21 Inverse the inverse of a square matrix D is denoted by D4 and it is the matrix such that D39ID I Vectors and Matrices Note that the vector 3 1 2 3 4 5 is really a 1 x 5 matrix We can use all of the matrix operations on vectors For example aT is a 5 x 1 matrix 5 rows 1 column or itis called a column vector If we have to vectors a 1 1 3 and b 2 4 5 then 2 abT 1 1 3 4 a b 21 14 3 5 17 5 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 2 Matrices 1 2 4 5 aTb 1 2 4 5 2 4 5 3 6 12 15 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 2 Matrices Matrix Operations with MATLAB In MATLAB matrices can be entered manually or by using some predefined MATLAB functions Here are some examples Examples A1 2 3 1 4 7 2 6 1 4 5 O This creates a 4 x 3 matrix each roWis separated by a semicolon entries Within a row are separated by a comma like vectors Bzeros4 This creates a 4 x 4 matrix of all zeroes Cones34 This creates a 3 x 4 matrix of all ones Deye5 This creates a 5 x 5 identity matrix Ediag123 This creates a 3 x 3 diagonal matrix With 1 2 and 3 on the first second and third rows respectively F5C scalar multiplication G3 4 5 1 6 3 2 2 9 1 O 17 HGA matrix addition JzGA matrix subtraction KGA multiplies the ij39h entry of each matrix LA transpose of A MGL matrix multiplication NLG matrix multiplication O1 2 3 1 6 O 6 2 14 Pinv0 inverse of O PO does it equal the identity matrix AdetP determinant of P Other MATLAB commands contour makes a contour plot of a 2dimensional array matrix bO550 makes a vector from O to 50 in increments of 5 O51050 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 2 Matrices MATLAB Exercise 2 21 Create a column vector in which the first element is 31 the elements decrease with increments of 4 and the last element is 9 22 Create the matrix shown on the board by using the vector notation for creating vectors with constant spacing when entering the rows ie do not enter the individual elements 23 Create the matrix A of the preceding problem then use colons to address a range of elements to create the following vectors a Create fourelement row vector named va that contains the second through fifth elements of the third row of A b Create a fourelement column vector name Vb that contains the elements of the sixth column of A 24 Give and demonstrate commands that a find the largest and smallest elements of matrix A and b add the elements of each row of A to yield a column vector of sums 25 Create a matrix B with diagonal equal to 2 and the first offdiagonal equal to 1 26 Sometimes we want to plot data that represents a function of two space dimensions x and y or of one space variable and time x and t We can plot the data as a contour plot isometric 3D plot or some kind of color or gray scale shading These routines usually assume that the data is available on a uniform orthogonal grid Let s plot some results for a regular grid Consider the function f sinx ywhereOSxS3ampOSyS4 Discretize it onto a uniform grid with spacing 01 in the x direction and 005 in the y direction 3181 2511 points This simulates the results of a numerical code with these dimensions spacings and size I purposefully picked a non square domain and different increments in the two directions in order to insure that you can handle these situations a Create the x and y vectors using linspace or the colon command b Using matrix multiplication of the x and y vectors create a matrix with the product of mi in the if11 entry c Calculate fsinxy using the matrix you just created 1 Contour the discrete results using the contour command Use a maximum value of 10 a minimum of 10 and a contour interval of 02 Your contour map should have the correct geometry 3 x 4 and aspect ratio Plotting the 3x4 domain on a square as is often done automatically by plotting programs is not acceptable Label the x and y axes and include a plot title e Make a mesh of the results using the mesh command Label the x y and z axes and include a plot title Use the correct aspect ratio for the x and y axes Make a surface plot of the results using the surf command Label the x y and z axes and include a plot title Use the correct aspect ratio for the x and y axes f New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 1 Vectors Vectors vector has magnitude and direction eg velocity specific discharge hydraulic gradient scalar has magnitude only eg porosity specific yield storage coefficient unit vector a unit vector is a vectorwith magnitude 1 Vectors are usually identified by a lower case letter To differentiate between scalar and vector variables vectors are identified in boldface such as a or with an arrow 13 above the vector name When writing by hand they are often indicated by Lmderline g or sometimes an overline E 39 some printed papersbooks also use these notations In mechanics i j k are unit vectors in the Cartesian x y and zdirections respectively Other notation for Cartesian unit vectors includes nx ny nz Example Suppose a is a 2 dimensional vector in particular let a 2 i 3 j Then a is a vector of magnitude 2 in the xdirection and 3 in the ydiIection In an xy plot it could be drawn as an arrow with its tail at the origin and its head at coordinate 23 But vectors aren t tied to an origin so a could also be drawn as an arrow with is tail at 55 and its head at 78 In other words the position of the tail is arbitrary but the magnitude and direction of the vector is fixed The vector 21 can also be written as 23 Vector Operations Define a01io2joz3k b Bl i 81 j 5 k where 01 02 03 81 6 and j are scalars Scalar multiplication c a 0 01 i c 0 j c 053 k where c is a scalar result is a vector Addition a b 01 81 i 02 82 03 g k result is a vector Subtraction a b 01 81 i 02 Blj 03 3 k result is a vector Magnitude l a l all all 032 12 result is a scalar DotProduct a39bal loz2 log 3lallblcos6 where 6 is the angle between the vectors result is a scalar ll New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 1 Vectors Arrays An array is an ordered list of data often referred to as a vector as we and Matlab will do Example 1 Suppose we measure water level in a well at eleven different times We can define a time array t that contains the 11 sampling times and a head array h that contains the 11 water levels So we might have t 0 5 10 15 20 25 30 35 40 45 50 in hours h 100 990 981 975 970 966 964 963 963 963 962 in feet Here we have two 1 ldimensional row vectors t and h Time is the independent variable and head in the dependent variable For each time in the t vector we have a corresponding water level in the h vector therefore the order of the data in each vector is important Arrays can be organized as rows as above or in columns If we take the transpose of t and h we get column vectors 0 100 05 990 t E 39 h E 45 963 50 962 Example 2 Suppose that we examine a transect across the landscape and measure topographic elevation z vegetation density d and soil moisture 6 at N locations ie N the number of locations Letx denote the horizontal position along the transect Then all three metrics vary with x that is they are fmetions of position x However we don t have complete functions only values at the N discrete locations We then represent this data set by four N dimensional arrays with corresponding data x x1 x2 x3 xm xN z 21 22 23 2m ZN d d1 d2 d3 dm dN e 191 62 93 gt61 51 This data could be organized in the form of a matrix with four rows representing the location and three metrics and with N columns representing the data X x1 x2 x3 quot 39 erl xN z 7 21 22 22 quot 39 Z N71 ZN d d1 d2 d3 d N71 d N e 61 62 63 quot 39 6N4 6N We describe matrix dimensions as the number of rows by the number of columns Here we have a 4gtltN matrix on the right hand side RHS and a 4gtlt1 column vector ofvectors on the LHS 12 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 1 Vectors BASIC MATLAB MATLAB comm and structure gt varname command Where gt denotes the command line prompt varname is the name of the variable to be assigned and command is the MATLAB command lfno varname is specified the default is ans lfthe semicolon is omitted the result is displayed Some commands contain calls to fmetions that require input variables For example the fmetion sqrt takes the square root of the number that is inputted to it All input variables are enclosed in parentheses order is important Input variable can be numbers or variables Comments can be added on separate lines or at the end of a line by preceeding it with the symbol Examples gt a 7239 gt ba5 l39 gt cb46al8abA239 gt dcosl 57sqrtcba39 gt e 123 this defines an array row vector e with elements 1 2 and 3 gt f 321 gt gef39 if two vectors have the same number of elements they can be added or subtracted using and 3 gth equot this transposes the row vector into a column vector Some Basic Commands Whos prints a list of all variables and their size help prints a list of help categories help command prints help for command eg gt help sqrt plot makes an xy plot of data requires input quit exits MATLAB size displays the size of the variable requires input linspace generates an array of evenly spaced numbers requires input title adds a graph title requires input xlabel adds a label to the xaxis of a plot requires input hold on freezes the plot so that additional data can be plotted on the same plot hold off unfreezes the plot comment line l3 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 1 Vectors MATLAB Exercise 1 a Build the vectors t and h in example 1 p 12 and plot h as a function of I Use a symbol for h Do not connect the symbols with lines Label the axes and title the graph b Build the vectors x z d and 0 in example 2 p 12 and plot 2 d and 0 as a fmetion ofx on the same graph Use different symbols for each dependent variable Connect with lines Label the axes and title the graph Note that example 2 does not assign actual numerical values we ll have to do that here Assume that the measurements are evenly spaced in x over a distance of 1km Use the matlab linspace command to make this vector Then use fmetions of x for the other parameters Make the elevation to be a linear fmetion of x so that z 100 01 xxo where x0 is the initial location and x is measured in meters Make vegetation density to be a sinusoidal fmetion of x with a wave length of 10m a mean value of 04 and an amplitude of 02 so that d 04 02 sin 27r10m Finally make moisture content a stochastic variable with mean value of 02 and a random component with a maximum range of 01 0 02 01 s The random variable is Lmiformly distributed between 1 and 1 and spatially uncorrelated Use the matlab rand command c You will often find occasion to examine the output of computer codes that you or someone else develops In some cases the data will consist of values of one or more variables as a fmetion of time I or space x With time and with onedim ensional spatial models the results will be presented as graphs of one or more dependent variables heads concentrations temperatures etc vs an independent variable usually time or distance The data will often be evenly spaced as a fmetion of the independent variable but that will not always be the case Data is usually in the form of a series of data pairs triplets etc consisting of a dependent variable value and its location say tor x This problem allows you to learn how to plot these types of graphs using MATLAB The structure is similar in other math codes and in plotting packages like PVWAVE Consider two fmetions of space x h cosx239 k x sinx39 Plot these two fmetions in the domain 0 S x S 3 using 31 points in x created with the linspace command Represent each fmetion by a different symbol Do not connect the points Provide a concise listing of the commands that you use to create the data file and plot it Label the x and y axes appropriately Include a plot title and legend a Use evenly spaced points in the interval 03 with Ax 01 increments b Place points spaced more closely at the larger x where there is more action using locations given by x with 0 g 139 g 9 and Ai 03 Compare to the plot in part a Note that both a and b use 31 points but b uses them more efficiently by employing uneven spacing of the independent variable Here you specified the region where more points will be used The most sophisticated plotting routines go one step further adaptively computing where more points are needed The symbolic math code Mathematica does this automatically 14 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming Writing Programs in MATLAB A MATLAB program is a file containing a list of MATLAB commands that are executed in a particular sequence In MATLAB a program file name must end with the extension called an M file To run a program type the file name without the m extension at the MATLAB prompt Each of the MATLAB exercises we ve done so far could be written in an Mfile The M file would contain the list of commands that you typed at the MATLAB prompt The advantages of using an M file instead of the comm and line are l the program can be run many times without reentering the commands 2 if you have errors it is easier to debug the program 3 the program is not removed from memory when you exit MATLAB and 4 control structures such as loops and if statements are easier to use inside of a program Programs should be specific enough to accomplish the required tasks but exible enough so they can be used more than once For example if you want to find the mean of five head values the program would be too specific if you coded the five values into the code A more flexible code would prompt the user to enter five head values and would calculate the mean of those values Then we could use the program to calculate the mean of any set of five head values In an even more exible code we could prompt the user for the number of head values that will be entered and then prompt the user to enter all of the head values This program could be used to calculate the mean of any set of an arbitrary number of head values Creating a MA TLAB Program To create a new MATLAB program select File New Mfile This will open the editor window Type your commands in this window When you save the file you should save it on your home directory which is U and save it with a extension To run the program either go to the command window and type the file name without the extension or select Tools Run from the menu bar in the editor window In either case MATLAB must be running in the directory where your file is located To find out what directory MATLAB is in from the command line type pwd at the prompt This is a L1an command to print the name of the current directory To change to the U drive type cd U to change directories type cd dimame where dimame is the directory name You can also edit an existing program by selecting File Open in the MATLAB menu bar Topdown Design Topdown design is a procedure for logically dividing the program tasks into the necessary MATLAB command statements At the top level we define the main tasks of the program Once these tasks are all defined we subdivide the tasks into smaller subtasks and continue this refinement Lmtil we are at the lowest level which is the MATLAB commands This provides a logical procedure from getting from the problem statement to the actual list of commands The problem statement might be very complex so we break it down into smaller subtasks that can be solved easily This procedure is very useful in large programs Similarity to other Programing Languages Many of the MATLAB programing conventions discussed below are very similar to those in the C and C languages and somewhat similar to those in FORTRAN For MATLAB it can be particularly useful to refer to a C or Ci r primers You can find links to these and to Unix primers at httpwww aintm arv d 39 39 html A particularly good C primer is at httpwwwcmcfacukDaveCCE html Keep in mind however that MATLAB works on vectors and matrices MATLAB is a vector language while C and C work on scalars 3l New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming FORTRAN is also a scalar language For vectorized supercomputers eg a old CRAY there are vector versions of these standard languages There are also special scalar versions meant for highly parallel computers such as an Intel Paragon Input and Output IO Input allows you to pass information into your program output allows you to get information out of your program To input data from the command window use the input command An example is Kinput Enter the hydraulic conductivity gt 39 The text inside the single quotes will be displayed in the comm and window and the program execution will pause Lmtil the user enters a value The value that is entered will be assigned to the variable Some variables are numbers either integers or real numbers other variables are character strings The command above is used to enter numbers Suppose we want to prompt the user to enter the name of a data file which is a character string We could use fnameinput Enter the name of the data file gt s 39 The s tells MATLAB that the variable zame is a character string Suppose we calculate the specific discharge q and want to display the results in the command window We have a choice of output commands including the following three q diSPOD fprintf The specific discharge is 82f q39 The first command is just a regular command statement without a semicolon so the results are displayed to the screen In this case the output will specify the variable name and its value The second disp command will display only the value to the screen39 it won t display the variable name The m39nzf command is used to display text along with the variable value In this case the text inside the single quotes will be displayed on the screen however the 82f will be replaced by the value of q The indicates the format of the variable whose value will be displayed using C language notation Some formats are f real number i integer e scientific notation s string The numbers between and f specify the total width 8 and the number of decimal places 2 for the display of the variable For example if q0l the displayed output would be The specific discharge is MM010 where the A indicates a blank space Other form at comm ands can be put inside the single quotes to put a tab in the output t or to start a new line after the output is displayed n For example fprintf The specific discharge is t 82f n q39 will tab between is and the data value and it will start a new line in the command window after the output is displaye d 32 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming Languages such as C and FORTRAN make extensive use of these kinds of formated input and output It is a good idea to become familiar with them here before taking up those languages For more on the format conventions used in MATLAB see a Mathlab primer or users manual or a C languarge primer and look for appropriate pages on formated output Programming Exercise Write a MATLAB Mfile that calculates specific discharge from hydraulic conductivity and hydraulic gradient Write the code so that it I prompts the user for the values of hydraulic conductivity and hydraulic gradient I uses these values to calculate specific discharge I displays the result Selection Structures Conditional iF39 statements Sometimes we want to execute a series of commands only if a certain condition is true In these situations we use an if statement The proper structure of an if statement is alt5 aa l 39 end This statement checks if the variable a has a value less than 5 If it does a is incremented by 139 otherwise nothing is done The condition that is tested is alt5 39 if this condition is true the commands after the if statement and before the end statement are executed Other relational operators are lt less than gt greater than lt less than or equal to gt greater than or equal to equal to N not equal to We can test two conditions simultaneously Suppose we want to increment a by 1 if it is greater than 5 and less than 10 We would use and represented by in the following statement if a gt 5 amp a lt 10 aal39 end For this condition to be true the value of a must be both greater than 5 and less than 10 Suppose we want to increment a by 1 if is it less than 5 or greater than 10 We would use an or represented by in the following statement ifalt5iagtlO aal39 end For this condition to be true the value of a must be either less than 5 or greater than 10 33 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming Sometimes we want to execute a group of commands if a condition is true and a different group of commands if the condition is false In this case we use an ifelse statement as follows if a gt 5 aal39 else aa 139 end In this example if the condition is true ie if the value of a is larger than 5 then a is incremented by 139 if the condition is false then a is decremented by 1 If we want to test multiple conditions we can use ifelseif statements as follows if a gt 5 aal39 elseif a gt 10 aa239 else aa 139 end For more on if statements and other selection statements below used in MATLAB see the Matlab primers or a users manual It can also be useful to refer to a C languarge primer Look appropriate pages on alternative selection statements Programming Exercise Write a MATLAB Mfile that calculates effective hydraulic conductivity of a twolayer system Write the code so that it I prompts the user for the values of hydraulic conductivity and thickness of the two layers I prompts the user to specify whether the arithmetic harmonic or geometric mean is to be used use a for arithmetic h for harmonic and g for geometric I calculates the appropriate conductivity using an ifelseif statement to determine which equation to use39 and I displays the results Run this program to calculate the effective hydraulic conductivity for K15 L1lO K2l L220 using arithmetic harmonic and geometric means where K is conductivity L is thickness and the subscripts indicate layer number Switch Statements Another selection structure is a switch statement Suppose a variable a can take on specific values of l 0 l and 239 and suppose that we want to execute different commands depending on the value of the variable We can use a switch statement to select among the several cases based on expression with a format as follows switc a case fl bsqrta239 dispb bsqrta239 34 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming dispb end Programming Exercise Repeat the last programming exercise but replace the ifelseifstatem ent with a switch statement Loops Sometimes we want to execute the same commands a known number of times For example we may want to write out the values of an nelement vector b If we know the size of the vector then we know the number of times we want to execute the commands We can use afar loop to repeat the commands a known number of times Suppose bl439239 then the number of elements is n5 We can use the followingfar loop for il n fprintf The value of b at time step i is i n ibi end In this statement we defined a loopcontrol variable 139 that takes the values of l234 and 5 in that order Since there are five possible values of 139 the loop is executed five times The output of this statement would be The value ofb at time step 1 is l The value ofb at time step 2 is 4 The value ofb at time step 3 is 3 The value ofb at time step 4 is 9 The value ofb at time step 5 is 2 The statement il n controls the number of times that the loop is executed and it controls the values of 139 during each execution lfwe used the statement il2345 instead we would get the same results lfwe used the statement inl l we would get the results printed in reverse order In afar statement it is even valid to have the loop control variable take on noninteger values however if we use the loop control variable to index an array as we did in the above example with bi we could run into problems because the array index must be a positive 1nteger Programming Exercise Write a MATLAB Mfile that calculates effective hydraulic conductivity of an nlayer system using the arithmetic mean You will need two nelement arraysone for hydraulic conductivity values and one for layer thicknesses Write the code so that it prompts the user for the number of layers initializes the hydraulic conductivity vector and layer thickness vector ie creates vectors of all zeroes and of the appropriate size uses afar loop to prompt the user for the hydraulic conductivity and layer thickness of each layer uses afar loop to calculate the effective hydraulic conductivity displays the result Run this program for three layers with K15 L15 K2l L2l0 K302 L315 While Statement Sometimes we want to execute the same commands for an unknown number of times Lmtil a certain condition is met For example we may want to generate random numbers between 0 and l Lmtil one of the random numbers 35 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming exceeds 08 We want to execute the same sequence of commands but we don t know how many times In this case we use a while statement as follows a039 na039 while a lt 08 arandl 39 nana 1 39 end In each execution of this while loop a random number is generated rand1 and assigned to the variable a This continues as long as the condition in the while statement is true ie as long as the cmrent value of a is less than 08 We count the number of executions of the loop by incrementing the variable ml by one each time the loop is executed Before the loop begins we define initialize a and mi because some value of a is needed just to test the condition to enter the while loop ie if a in not defined we can t check if the condition alt08 is true We must be sure that the condition in the while statement will eventually be false otherwise the while loop will never stop executing and we will be in an infinite loop Programming Exercise Modify the previous programming exercise that calculates effective hydraulic conductivity of an nlayer system Instead of prompting the user for the total number of layers ask the user if he wants to enter information for another layer Then continue prompting the user for hydraulic conductivity and thickness information Lmtil the uses indicates he is done Then calculate the effective hydraulic conductivity using the arithmetic mean For this problem you can t set the arrays to the appropriate size up front because you don t know their lengths But you can initialize them to empty K39 Although this may not be necessary it eliminates any possibility of reusing leftover data from a previous run Array Indexing When using indexed arrays in for loops make sure that the arrays are indexed properly A onedimensional array ql 12 14 16 18 20 has six elements with q1139 q21239 q31439 q41639 q51839 and q620 so qn picks out the n3911 element of q If we had afar loop running from jl 7 and inside of the loop we used qj we would get an error index exceeds matrix dimensions and the program would stop running Note that if we had a for loop running from jl 22 thenj would take the values of 1 12 14 16 18 and 20 Ifwe used qj inside of the loop we would get an error because q12 does not exist That is there is a value of 12 in the q vector but there is not 12nd value of the q vector The array indices must be positive numbers PROGRAMING EXERCISES 1 Write a MATLAB Mfile that calculates specific discharge from hydraulic conductivity and hydraulic gradient Write the code so that it I prompts the user for the values of hydraulic conductivity and hydraulic gradient I uses these values to calculate specific discharge displays the result 2 Write a MATLAB Mfile that calculates effective hydraulic conductivity of a twolayer system Write the code so that it prompts the user for the values of hydraulic conductivity and thickness of the two layers prompts the user to specify whether the arithmetic harmonic or geometric mean is to be used use a for arithmetic h for harmonic and g for geometric calculates the appropriate conductivity using an ifelseif statement to determine which equation to use39 and displays the results 36 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Lab Notes 3 Programming 3 Repeat the last programming exercise but replace the ifelseifstatem ent with a switch statement 4 Write a MATLAB Mfile that calculates effective hydraulic conductivity of an nlayer system using the arithmetic mean You will need two nelement arraysone for hydraulic conductivity values and one for layer thicknesses Write the code so that it prompts the user for the number of layers initializes the hydraulic conductivity vector and layer thickness vector ie creates vectors of all zeroes and of the appropriate size uses afar loop to prompt the user for the hydraulic conductivity and layer thickness of each layer uses afar loop to calculate the effective hydraulic conductivity displays the result Run this program for three layers with K15 L15 K2l L2l0 K302 L315 5 Modify the previous programming exercise that calculates effective hydraulic conductivity of an nlayer system Instead of prompting the user for the total number of layers ask the user if he wants to enter information for another layer Then continue prompting the user for hydraulic conductivity and thickness information Lmtil the uses indicates he is done Then calculate the effective hydraulic conductivity using the arithmetic mean For this problem you can t set the arrays to the appropriate size up front because you don t know their lengths But you can initialize them to empty K39 Although this may not be necessary it eliminates any possibility of reusing leftover data from a previous run 37 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Note that for the step change problem T0t 05 To for tgt 0 The step smears over time and unlike the diffusion problem the concentration at the origin changes It is not a boundary condition Transform Solutions to Heat Diffusion see Crank 22 It can be shown by dimensional analysis that solutions to 30 are often of the form A x2 Txt rm 41 68 where I is a function to be determined In particular we ve already found Gaussian and error function solutions Notice that the function is self similar That is the function is stretched one way or another but the shape remains unchanged over timespace You ve seen such similarity solutions before in the Theis Well function which has similarity variable ur S4Tt where the parameters take on the normal wellhydraulics definitions and hydraulic diffusivity is DTS It is not unusual to assume a solution of the form of 68 and then to seek the coefficientl that matches the boundary and initial conditions As a demonstration let s take a look at one such application of similarity to find another path to a familiar solution Boltzmann Transformation Crank 1975 p 105 this page is not included in the supplemental reading Reconsider the diffusion equation given by 6T 6 6T D T 69 6t 6x 6x with initial and boundary conditions for a semiin nite domain 0 S x S 00 described by IC Tx0 T a constant BCl T0 t To a constant 79abc BC2 Too t T a constant From our previous work we might expect that the solution involves an error function Instead let s assume only that it is a similarity solution of the type shown in 68 De ne a similarity variable w 71 Then apply this change in variables and some chain rules to convert 69 into an ODE Converting a PDE to an ODE is the goal of many methods that solve parabolic PDEs In this case the chain rules are 19 Typically its some combination of exponential or error functions l99 New Mexico Tech Hyd5 Hydrology Program Quantitative Methods in Hydrology zav wzyImw x m at x at dn 2DlZt32 dn 4D12t32 dn 61 w 1 w 6x2 6x2 6177 I4Dt 6177 620 161 116nm111md20 6x2 2 6x 6x2 6177 t 14Dt 6177 6x2 6177 1 4Dt d7 14Dt d7 4Dt 61772 72abc Notice the ordinary derivatives with respect to the similarity variable 17 We can substitute into the diffusion equation aT x dT 62T DdzTor 6t 4D12t3 d7 6x2 2 4Dt any x 11312 4D12t32 d7 4t 61772 This simpli es to 2 2 d 2x d T0 or d 2nd T0 73 d7 J4Dt 6177 d7 6177 which is a second order ODE We can solve by reduction in order say be de ning dT v E 75 where the symbol v was chosen arbitrarily The ODE becomes 1st order i 2 77 v 0 76 6177 Its solution by separation of variables SOV and integration is l611v 2776177 gt lnv 772 c1 gt V a1 eXp772 v d T a1 exp 772 77 6177 where cl and al are constants and we have derived a new 1st order ODE to solve We can apply SOV again to get 200 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology dT a1 exp 77sz 78 Obviously we should integrate this expression However to progress further we need to apply initial and boundary conditions Recall that it is these conditions that determine the form of the solution to the PDE In fact with the Bolzmann approach it is only when the initial and boundary conditions are expressible in terms of I alone and x and t are not involved separately that the Boltzmann transformation can be used Crank 1975 p 106 For example in our case when t0 no matter the value of x then 1700 Similarly when x0 then 170 What about when xoo We have 7oo again Thus our conditions become IC amp BC2 T7oo T a constant BCl TI0 To a constant 79ab Let s use the second condition and integrate from 170 to 17 4Dt T x LU dT alLmexp 77sz T To a1 Eerf 43 80 where we can recognize the integral as appearing in the error function That leaves one constant to determine from the other condition T07 00 T where erf i7 erf oo 1 or 7239 2 T To alger w a a 70 To 81 Substituting back into 80 gives the nal solution of the PDE T T0 erf x 82a T T0 4Dt Flux calculation at x0 One application of this solution is to the calculation of the time varying ux of heat into the domain from the boundary or to the calculation of the total heat added The rst of these involves Fourier s law and taking the derivative of 82 wrt x while the second involves the integral of that derivative over time or of 82 over space The flux is given by x TxtTZ T er 82b 0 0 fm dT d x D39ff39 K KT T 83 1 us1ve ux dx 1 0dx erf where K is the heat conductivity and DKcp The derivative on the right evaluates using the chain rule p 24 of these notes and by applying the derivative of an integral with a variable upper limit ie the derivative of an error function p 44 of these notes 201 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology dT i 1 dx T To dx er 07 T To dx d77erf77 d x 2 2 T To E m leP U 84 To x2 t eXp41 Thus the diffusive ux any xt is Z Flux K T Kwexp x 85 x V 7th 4Dt At the boundary the flux into or out of the domain is given by applying this model at x0 0r dT T T T T Bounda d1ffus1ve ux at x0 K K K c 86 W dx Dt J p m 160 The boundary flux is greater for higher thermal conductivity and heat capacity c larger temperature differences T0 Z and smaller time It decreases with the square root of time We integrate this over time to get the quantity of heat added as a function of time H t or To T 39 z MLIzzf H JKpc J 1KPC J Int KPCT0 TIMZ 87 The amount of heat added increases with the square root of time The amountH should also equal the spatial integral of heat added or H prchxr T11dx pom 4 1 era16k pom T flerfc x Me pom Tux41 I erfc77d77 88 V4Dt wKpc To T 7139 which agrees with 87 I had to look up the integral at 39n n quot id wolfram com Frfc html the mathematica web site to nd that I erfc77d77 f 89 In application if I had not been able to nd this integral I would have relied on 87 It s nice to have more than one way to solve a problem 202 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Other applications We can apply this model elsewhere in hydrology One common application is bank storage due to a passing ood wave If his the initial head in the aquifer ho is the height of the ood wave using the same datum D is hydraulic diffusivity TSy then by homology to 82 86 and 87 the aquifer head the boundary ux into the aquifer and the volume of water stored are given by x hxt h hl h erf J4Dt dh h h h h T T 0 TS 90abc Fm V y 5 VJTSyh0 h1 E 7239 This tells us that the rate of bank infiltration and the amount of bank storage is greater for higher transmissivity and storage coef cient and larger head differences h0 h The rate decreases with the square root of time while the storage volume increases with the square root of time Another common application is moisture absorption into a drier or desorption from a wetter porous material due to capillarity and ignoring gravity provided we assume that capillary diffusion is a linear process If 9 is the initial moisture content in the soil volume of water per volume of soil 90 is boundary moisture content D is capillary soil moisture diffusivity K dlldlg where 1p is soil moisture tension and K is hydraulic conductivity then by homology to 82 86 and 87 the moisture content boundary ux into the soil and the total amount of absorption are given by 6xt 60 02 90 erf DJ 91abc 43 dx x0 7113 4 V 5090 a f 7139 The rate and amount of moisture absorption is greater for higher soil moisture diffusivity and larger moisture differences 60 Q The rate decreases and the amount increases with the square root of time If we add gravity20 to this moisture problem we can study in ltration but continuing to assume that K and D are constants and the problem linear In this case only the flux changes or DKXOJDK 92 Thus the presence of gravity increases the in ltration rate but the rate is suppressed by higher initial moisture 19 The actual problem is more complicated than this Both D and K depend on 20 Assuming thatx is positive downward 203 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology moisture content The problem is actually nonlinear Let s return to the heat diffusion problem to explore this issue Nonlinear diffusion Reconsider the diffusion equation but for the case of a diffusion coefficient that depends on the state For heat diffusion this would be given by 6T 6 6T 5 ing 93 Other applications of state dependent diffusion coef cient include solute diffusion groundwater hydraulic diffusion where in a phreatic aquifer the transmissivity depends on the saturated thickness and thus on the water table elevation and moisture absorption and in ltration How can we solve 68 for any general function DT without specifying that function explicitly We ll see that what we can do is to develop a solution approach implement it and carry out the symbolic calculations to a certain point but then we ll have to apply a particular function DT and complete the remainder of the calculation numerically Below we assume thatD depends on temperature through a dependence of thermal conductivity KT on temperature Heat capacity is assumed constant While there are some very good nonlinear solution approaches21 for 68 also capable of handling advection we ll focus on the simplest It is to illustrate only you won t be applying this method yourself Basically the idea is to use the Boltzmann transformation again We assume that the solution will be something like that in the similarity solution of 68 apply the change of variables in 71 but in terms of a fixed D0DT0 taken at the boundary x 94 1 4D0t and derive a nonlinear ODE similar to 73 or 1DTdT 2 dT0 95 6177 D0 d7 6177 This is presented in some detail in Crank 1975 Chap 7 and reviewed in the context of soil moisture absorption and infiltration in Eagleson Dynamic Hydrology Prentice Hall 1970 pp 291 299 For example the boundary diffusive ux at x0 becomes 1 KPCLT x0 2 J a777 21 John Phillip CSIRO and Yves Parlange Cornell are well known in hydrology soil physics and mathematical physics for their contributions along these lines For example both worked on something called the flux concentration approach KdT dx 86 710 204 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology which by comparison to the linear case we can rewrite as KdT dx H Jpr 87 where the weighted mean conductivity is given by 2 z 1K0 2 d1 4 T To dn and K 0 is taken at the boundary K0KT0 chT0 Crank Conduction of H eat in Solids Oxford Press 1959 did some numerical experiments to determine approximations for I and then explored the solution for various functions KT You can also nd some of this in Eagleson 1970 For example Crank found that as long as K and D increase with Tthen the best weighting for increasing temperature cases is given by I 2 gm T WI T Z WWW lt89 88 710 One would take experimental data on K T perform this integration perhaps numerically and insert into 87 to get the model for heat ux Similar calculations are done for temperature distribution and total content Homologous calculations are done for stream bank storage and moisture absorption and in ltration 205 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology SecondOrder Linear ODEs Textbook Chap 2 Motivation Recall from notes pp 5859 the second example of a DE that we introduced there d2 abc2 1 1 Zt ZQFQW 31 This equation represents conservation of water mass actually conservation of water volume the model assumes an essentially constant water density in a con ned but leaky aquifer The unknown is denoted by the letter L representing hydraulic head and the independent variable by the letter x L representing location It is a second order ordinary differential equation ODE The symbols A 0 and Qw denote known parameters Respectively they represent aquifer transmissivity divided by a leakage coefficient l L2 TB K see p 59 and the paragraph on dimensional homogeneity the head L in a vertically adjacent overlying or underlying aquifer and the pumping per unit area Lt All terms in a have units of lL To avoid confusion with the notation in the text which uses the symbol 1 for something else let s replace l with the symbol l or 1 Qw39 a2 Another DE presented on p 46 represents solute mass balance in a owing river with the time rate of change of mass in storage balanced by advective and dispersive uxes or 0 0 Here the unknown is denoted by the letter C representing solute concentration There are two independent variables denoted by the letters x space and t time As there are now two or more independent variables it is a partial differential equation PDE The highest order derivative in x defines it as a second order equation in space while the highest order derivative in t de nes it as a rst order equation in time The symbols v and D denote parameters the mean river velocity and a dispersion coef cient Suppose that there is a continuous unchanging source of solute and that the ow rate in the river is constant Then aCat0 the solute concentration is in steadystate and b becomes a second order ODE in x 2 Dd C dx dxz 0 C The rst term represents solute advection and the second dispersion which are in balance Notice that the partial differences 6 have become ordinary differences d s as there is now 120 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology only one independent variable x Suppose further that the solute is organic and that it can biotransformed by in situ reactions The simplest model of this is that the solute concentration decays exponentially with time which in c is modeled by adding a socalled rstorder decay process 2 vZ C DZX aC d1 x where a 1 T is the decay constant Here advection and dispersion are balance by decay We can write this in mathematical standard form by dividing by the quantity 7D and gathering all terms on the lefthand side 2 d L C0 d2 dx D dx D This is secondorder ODE in spatial location x These are two common examples of secondorder ODEs encountered in hydrology but there are many others These equations are also linear if their parameters are independent of the dependent variable They can vary with location x however Introduction Textbook Section 21 These two examples belong to the general class of linear secondorder ODEs which your book writes as 1 y pxy qxy rx where r is the forcing and for a linear model parameters p q and r do NOT depend on the dependent variable y but can depend on the independent variable x How do the two examples above t this model Model Translation to egn 11 51 1 1 1 1 E oQw39 y axxaP05qEqarE oQw C0 yCxxp 1q a 2 r0 dx Ddx D D D Note that p q and r can vary with x and the models remain linear Thus 0 and Qw in the rst model and v D and 0c in the second model can also vary in x for linear models New Mexico Tech Hyd 5 Hydrology Progmm Quantitative Methods in Hydrology Then from the text Chapter 2 If rtx E 0 that is Mr 0 for all X considered read r39 is identically zeroquot then 1 reduces to 7 rquot 12thquot 1m 0 and is called homogeneous If r SE 0 then l is called nonhomogeneous This is similar to Sec 15 For instance a nonhomogeneous linear ODE is y 6quot cos x U n y 2 and a homogeneous linear ODE is l n39 iquot x v O in standard form rquot i l y 0 An example of a nonlinear ODE is y quot2 0 The functions 1 and q in l and 2 are called the coef cients of the ODEs Solutions are defined similarly as for first order ODEs in Chap 1 A function x an is called a volniion of a linear or nonlinear secondorder ODE on some open interval 1 if 11 is defined and twice differentiable throughout that interval and is such that the ODE becomes an identity if we replace the unknown y by Ii the derivative quot by I7 and the second derivative 7 by II Examples are given below Fundamental Theorem for the Homogeneous Linear ODE 2 For a llulllagcnmax linaai ODE 2 any linear owninalinn affirm 39olntiun39 In an open inmn39al 1 ix again a wlnrinn gill an I In particular for II39I an equation nuns and mnxianf multiples ifmlntianx are again mniinns Not in the text s Chapter 2 For a secondorder homogeneous linear ODE 2 a boundary value problem consists of equation 2 taken over a finite interval a S x S b with two boundary conditions BCs one at each end of the interval The boundary conditions can be written as k1 3 01 2 01 Ka 11 yb lzy39UJ Kb 6 122 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology where the k s and 5 are parameters weighing the relative contribution of y and y at each boundary and the K s are values For example if you know the value of y at a then In 1 kg 0 and Ka is the value ya In PDEs this is called a first type or Dirichlet boundary If you know the gradient of y at b then then 11 0 12 1 and Kb is the value y b In PDEs this is called a second type or Neumann boundary Finally at a if both In and k2 are nonzero it is called a third type or Cauchy boundary in PDEs All three boundary condition types are used to solve equations of the type given in 1 However Chapter 2 of your textbook focus es on initial value problems for 2nd order ODEs In hydrology these problems are less common than 2nd order BVPs but they are very common in other areas of mathematical physics like geophysics Back to the text For a second order homogeneous linear ODE 2 an initial value problem consists of 2 and two initial conditions 4 390 K0 quot 0 K1 These conditions prescribe given values K0 and K1 ol the solution and its rst derivative the slope 01 its eun39e at the same given 39 390 in the open interval considered The conditions 4are used to determine the two arbitrary constants 391 and 392 in a general solution 5 J 1quot1 quot2quot2 of the ODE here 1 and 2 are suitable solutions of the ODE General SolutionY Basis Particular Solution A general solution of an ODE 2 011 an open interval is a solution 5 in which 1 and 72 are solutions of 2 on I that are not proportional and 1 and 2 are arbitrary constants These 1 72 are called a basis or a fundamental system of solutions of 2 on I A particular solution of 2 on I is obtained if we assign specific values to 1 and 2 in 5 Or for the BVP the conditions in eqn e on the previous page are used to determine the two constants c1 and CZ Two conditions are needed as this is a 2 01 order equation In the IVPs these two initial conditions 4 are taken at the beginning of the domain and in BVPs these two boundary conditions e are taken at each end of the domain In either event the general solution is given by 5 123 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Two l unctions 1 and 392 are called linearly independent on an interval I where they are defined if 7 kLi39ltx k2y239 O eveiywhere on I implies k1 and k2 0 And 1 and 2 are called linearly dependent on I if 7 also holds for some constants k1 k2 not both zero Then it k1 0 or kg at t we can divide and see that 1 and 2 are proportional quotr r or T k li 1 i 2 2 i 1 kl k2 In contrast in the case of linear fill INHIL IFIICI these functions are not proportional because then we cannot divide in 7 This gives the following Basis Reformulated A basis 01 solutions 01 2 on an open interval I is a pair 01 linearly independent solutions of 2 on I Reduction of Order Section 21 Reduction of Order if a Solution ls Known Basis Find a basis ot solutions 01 the ODE x2 7 Myquot 7 taquot Ar 0 Solution Inspection shows that i39l x is a solution because j39l and r 0 so that the first term vanishes identically and the second and third terms cancel The idea of the method is to substitute 1 r n u i y Iri39l in l39 I i39 u y n t39 211 into the ODE This giVes 2 7 39Iquotr 24 7 rtu x HI M39 0 Mx and 7in cancel and we are left with the following ODE which we divide by x order and simplify r2 7 XVIIIX ZIIII 7 2II 0 39Z 7 luH r 7 zm 0 This ODE is ott irsl order in t H namcl t2 39ltrI t ll tr 2 0 Separation of variables and integration giVes 1 2 7 7 7 tl In Ivl In I 7 II 7 2 In In 7 l We need no constant ol integration because we want to obtain a particular solution similarly in the next integration Taking exponents and integrating again we obtain I 7 z 7 7 II It IN In T hence 392 II 39 In I Since and r In I are linearly independent their quotient is not constantL we have obtained a basis of solution valid for all posiln c 124 New Mexico Tech Hyd 5 Hydrology Progmm Quantitative Methods in Hydrology SecondOrder Linear ODEs with Constant Coefficients Section 21 We shall now consider secondorder homogeneous linear ODEs whose coefficients a and b are constant 1 yquot aquot by 0 These equations have important applications especially in connection with mechanical and electrical Vibrations as we shall see in Secs 24 28 and 29 How to solve l We remember from Sec 15 that the solution of the firstorder linear ODE with a constant coefficient k quot k O is an exponential function y 60 This gives us the idea to try as a solution of l the function Substituting 2 and its derivatives 7 M and y y into our equation 1 we obtain A2 mt by Hence if L is a solution of the important characteristic equation or auxiliary equation 3 A2ztb0 then the exponential function 2 is a solution of the ODE 1 Now from elementary algebra we recall that the roots of this quadratic equation 3 are 4 A1 7a V112 e 4b A2 w e V512 7 4b 3 and 4 will be basic because our derivation shows that the functions 5 1 8A1 and 2 2 are solutions of 1 Verify this by substituting 5 into 1 From algebra we further know that the quadratic equation 3 may have three kinds of roots depending on the Sign of the discriminant 2 7 4b namely Case 1 Two real room ifa2 7 411 gt 0 Case 11 A real double r001 ft2 7 417 0 Case 111 Complex conjugate roots If a2 7 41 lt 0 125 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Two distinct real roots Case I Two Distinct Real Roots A and A2 In this case a basis of solutions of l on any interval is h A1 and 2 9 because 71 and 72 are defined and real for all X and their quotient is not constant The corresponding general solution is 6 7 la 1 626 Initial Value Problem in the Case of Distinct Real Roots 0ch the initial value problem 7 y 7 239 0 n0 4 MO 754 SOIlltiUlI Step I General snluriou The characteristic equalion is A2 A 2 0 115 roots are A1l l 31 and A2 l 5 2 gt0 thal we obtain the general solution r z lv39T 20 21 126 NeWMexico Tech Hyd510 Hydrology Program Quantitative Methods in Hydrology Step 2 Particular solution Since quot39l rip lt39zCAQI we obtain from the general solution and the initial conclilionx Vl39lOl 1 392 4 V39l0l 1 lt39z 5 Hence 1 2 and 392 3 This ghes the IIxn39z39ry quotr 38 h Figure 29 lthon39s that the curve heginx at V39 4 with a negative slope 175 but note that the axex hme different walesli in agreement with the initial conditions 0100 ON 0 05 1 15 2 139 Fig 29 Solution in Example 2 Real double root Case II Real Double Root A a2 If the discriminant a2 7 4b is zero we see directly from 4 that we get only one root t tl A2 12 hence only one solution Y1 Beta2M Hence in the case of a double root of 3 a basis of solutions of l on any interval is Fax2 Kean2 The corresponding general solution is 7 N 1 2 c39 quotquot2 127 New Mexico Tech Hyd 510 Hydrology Progmm Quantitative Methods in Hydrology Two complex roots 1 1 Case lll Complex Roots a nu and 3a 1 This case occurs if the discriminant 2 e 4b of the characteristic equation 3 is negative In this case the roots of 3 and thus the solutions of the ODE 1 come at first out complex However we show that from them we can obtain a basis of real solutions 8 1 fax2 cos wx 2 e fax2 sin tux w gt 0 where 42 e in It can be verified by substitution that these are solutions in the present case We shall derive them systematically after the two examples by using the complex exponential function They form a basis on any interval since their quotient cot an is not constant Hence a real general solution in Case 11 is 9 y Fax2 A cos aux B sin an A B arbitrary Summary 2nd order ODE with constant coef cients y ay by 0 Characteristic eqn 23 a1 b1 0 Summary of Cases llIl Case Roots of 2 Basis 0f 1 General Solution of l Distinct real I 1quot quot21 39 KIM 39zz AZ Al A2 39 Real double root N I A 71 ram2 rye 1112 y 1 Ewe L312 7 i Complex conjugate fax2 cm a 111 A1 a in Wl2 quot 39 far2m cos u39 B sin arr A 1 I 39 sin cox 2 u Iw Z a where ml bif andA B cl and 02 are constants 128 New Mexico Tech Hyd 5 Hydrology Progmm Quantitative Methods in Hydrology Nonhomogeneous ODEs Section 27 Method of Undetermined Coefficients In this section we proceed from homogeneous t0 norihnmogeneous linear ODEs l 7quot 1xv qxv rx where rx SE 0 We shall see that a general solution of l is the sum of a general solution of the corresponding homogeneous ODE 2 mm aim1v 0 and a particular solution of 1 These two new terms general solution of l and particular solution of 1quot are defined as follows General Solution Particular Solution A general solution of he nonhoniogcncous ODE l on an open interval 1 is a solLIlion ol the form 3 TU MAX TAX here I 2 Ln avg is a general soluliou of the homogeneous ODE 2 on and p is any solution of l on containing no arbitrary constants A particular solution of l on I is a solution obtained from 3 by assigning specific values to the arbitrary constants 391 and 2 in y To solve the nonhomogeneous ODE 1 or an IVP or BVP for 1 we have to solve the homogeneous ODE 2 for a solution and then find m solution yp of 1 to obtain a general solution 3 of equation 1 See the book for details and proofs of why this works The method of undetermined coef cients There are more general methods to find yp but we will focus on a simple method that works for most problems that you are likely to encounter The more general method is Variation of Parameters and is discussed in 220 not covered in lecture The simpler but restricted method is called the Method of Undetermined Coef cients The basic idea is to guess all possible terms in the solution yp then solve for their coefficients term by term to satisfy 1 The result is yp The guess consists of sums of terms like those in rx This works only if rx is composed of terms involving simple functions in particular powers of x trigonometric functions and exponentials or their products 129 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology More precisely the method of undetermined coefficients is suitable for linear ODEs with constant coefficients a and b 4 y ay by rx when rU39 is an exponential function a power of x a cosine or sine or sums or products of such functions These functions have derivatives similar to I X itself This gives the idea We choose a form for yp similar to rx but with unknown coefficients to be determined by substituting that yp and its derivatives into the ODE Table 21 on p 80 shows the choice of yp for practically important forms of rx Table 11 Method of Undetermined Coefficients Term in its Choice for yp39 ATV 4 quotCYV hitquot it 0 KM l1 1 Kls K0 it cos w39 JKcos or M sin or lt Sm LOX he cos a39 I M e quot l K cos at M sin tux e 39 Slll39 Choice Rules for the Method of Undetermined Coefficients a Basic Rule fr in 4 is one of the functions in thefirst column in Table 21 choose yp in the same line and determine its undetermined coef cients by substituting yp and its derivatives into 4 Modi cation Rule If a term in your choice for yp happens to be a solution oft1e homogeneous ODE correspomling to 4 multiply your choice ofyp by X or by 1392 if this solution corresponds to a double root of the characteristic equation of the homogeneous ODE c Sum Rule Ifrx is a sum offunctions in the first column of Table 21 choose for yp the sum of the functions in the corresponding lines of the second column A 739 v The Basic Rule applies when rx is a single term The Modification Rule helps in the indicated case and to recognize such a case we have to solve the homogeneous ODE first The Sum Rule follows by noting that the sum of two solutions of l with r r1 and r r2 and the same left side is a solution of l with r r1 r2 Verify Your book appears to neglect products but you can handle products of the basic functions in the rst column of Table 21 You just apply the modi cation and sum rules to the product 130 New Mexico Tech Hydrology Progmm Hyd 510 Quantitative Methods in Hydrology Application of the Modification Rule b Solve the initial value problem 0 y 3H 225y 7 7 I0 7151 0 7 I M0 7 0 Solution Ste 1 General solution 0 the llama eneous ODE The characterix39tic etuatioii of the P g l lltHnOgCllCOUh ODE is A2 3A i l 2 0 Hence the homogeneous ODE has the general solution yh 1 t392i39oT1395T Step 2 Solution y of the Imnhamageneaus ODE The function TI39ST on the right would normally require 7 r the choice 7 1quot But we see from An that this tuiiction IS a SOlllllOl ot the homogeneous ODE WlllCll corresponds to a dimIto mm of the characteristic equation Hence according to the Modi cation Rule we have to multiply otir choice function by X2 That is we cho se 028715 1quot Then C2 7 I t Zie lSl39 yg 7 2 7 3x 7 339 2 iizwlSx We subatitute thelte expressions into the given ODE and omit the actor 715139 Thigt yields 2 7 6 22512 3cm 7 15 22502 7 7 i0 coefficients ol39 390 39J39l Hence the given Comparing the give 0 0 0 0 2C 710 hence C 75 This gives the solution 2 71 39 e yp 73A DE ha the general solution Iquot yh Ai39p cl clezfl39sx Step 3 Svlil uli of the initial l39llllt 1l39lllt lll Setting t in 39 tlltl thing the I ii39xl initial Cl llLllllUlL we obtain V39l 391 I Differentintion of 39 ghei y 7 2 7 Sq 7 152 i515 7 iota 15139 75 21715quot r From llllm und the wound initial L Olltlllan we have y ill c2 7 l W1 0 Hence 1392 libt39l 15 This givelt the anwwer Fi 7 I 15 42 15 7 5X28 1395x 7 t 15x 7 SAM 1395 The curve beginx with it horizonliil tiii ont cranes the Hixia nt 0627 llL l39L l5 7 5 U and npproachex the zixilt from below as 39 inc 13gt I y 10 05 0 i 1 705 710 Fig 50 Solution in Example 2 131 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology Equidimensional ODEs Section 25 Euler Cauchy equations4 are ODEs of the form 1 A2 rmquot by O with given constants I and b and unknown 3939 We substitute 2 i39 quotquot and its derivatives 3 quotMW 1 and yquot HUN 7 lquotquot2 into 1 This gives 392mm 7 lr 2 Lt39mxm l bx t We now see that 2 was a rather natural choice because we have obtained a common factor X quot Dropping it we have the auxiliary equation mtm 7 l in b U or 3 III2 u I I 0 Note I 1110111 Hence 39 quotquot is a solution of 1 if and only if m is a root of 3 The roots of 3 are Case I If the roots ml and quotI2 are real and different then solutions are Mix quot and 2 sz They are linearly independent since their quotient is not constant Hence they constitute a basis of solutions of l for all 39 for which they are real The corresponding general solution for all these i is 5 y 1l39m 2 m392 r1 2 arbitrary Case II Equation 4 shows that the auxiliary equation 3 has a double root ml l 7 a if and only if 1 7 12 7 4b 0 The EulerCauchy equation 1 then has the form 6 iZy Mquot in A 50 0 Thus in this quotcritical casequot a basis of solutions for aosilive r is l39 r quot and 1 quot2 quot In x where m it 7 1 Linear independence follows from the fact that the quotient of these solutions is not constant Hence for all r for which 1 and quot2 are defined and real 1 general solution is 7 y 1 c2 ln3939quot m t a Case III The case of complex roots is of minor practical importance and it suffices to present an example that explains the derivation of real solutions from complex ones 132 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Series Solutions Chap 5 So far we ve examined the general class of linear secondorder ODEs which your book writes as 1 y pxy qxy rx where r is the forcing and parameters p q and r do NOT depend on the dependent variable y but can depend on the independent variable x The methods we ve examined have considered relatively simple and special functions for the parameters but what happens if the parameters have more complicated form If rx is more complicated than can be handled by the Method of Undetermined Coef cients 27 then we typically resort to the Method of Variation of Parameters 2 10 which we have not covered in glass But there is another method that works very well in this case called LaPlace Transforms Chap 6 which we use on 1st and 2quotd order linear equations and which is very versatile It even gives the particular solution right away without solving for the homogeneous solutions rst as an intermediate step We ll come back to it later when we study PDEs There is another method we use if px and qx have more complex functional form and its particularly valuable for BVPs In this approach we resort to series solutions We know that we can approximate functions by series The idea is to assume that the solution is written as a power series to plug that series into the equation and to equate sums of coef cients of equal power to zero solving for the coef cients in the power series For a second order equation each coef cient can be written terms of another coef cient resulting in a recurrence relation for the coef cients and leaving only two underdetermined coef cients to be resolved by boundary or initial conditions Often the series solutions are represented by special functions with new functional names such as the Bessel Function of first type of order zero Certain forms of 2quotd order ODEs occur often in mathematical physics engineering and hydrology The series solutions for these problems have been worked out in detail and are typically represented by the special functions It is useful to learn to recognize them or at least to know that you should look to see if the equation you are dealing with is of one of these types For example if you have a problem in radial ow such as in groundwater well hydraulics Bessel Equations are common while Legendre Equations are common in spherical ow or solute diffusion within a particle like a spherical microporous feldspar grain Later we ll see that the methods can be combined to solve PDEs For example consider a transient groundwater ow problem that is onedimension in space It can be LaPlace Transformed in time to remove time and reduce it to an ODE in space which is solved by one of the methods we have been studying An example is radial ow to a pumping well in an in nite leaky con ned aquifer After the LaPlace Transform in time the resulting ODE in LaPlace Space is a Bessel Equation for which one can look up and apply the solution 133 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology The series solution approach 39 39 39 from the text p 168 For a given ODE 1 pm thu 0 we first represent pr and x by power series in powers of x or of X 7 X0 if solutions in powers of x 7 0 are wanted Often I39 and q39 are polynomials and then nothing needs to be done in this first step Next we assume a solution in the form of a power series with unknown coefficients 3 3 y E amx 10 11x IZXZ 13x 3 14170 and insert this series and the series obtained by termwise differentiation cc a y E mamxnk1 11 My 313392 m7 4 x b yquot E mm 7 lumxm 2 212 3 2 13x 4 3I4XZ quot172 into the ODE Then we collect like powers of r and equate the sum of the coefficients of each occurring power of x39 to zero starting with the constant terms then taking the terms containing l39 then the terms in A2 and so on This gives equations from which we can determine the unknown coefficients of 3 successively The text then goes on in 52 to talk about theory If the series converges everywhere the function is said to be analytic If it diverges in certain locations those are called singularities A Real Analytic Function A real function ftx is called analytic 1 I point 39 I NO if it can be represented by a power series in powers of 39 0 with radius of convergence R gt 0 simple example is a divide by zero The function and its derivatives must converge to be analytic Convergence is de ned in the usual ways e g by comparing successive terms in the series Series solutions for analytical functions are straightforward it is all the special cases involving singularities that make solving a problem dif cult For certain of these more complicated equations the Method of F robenius 54 is particularly useful In any event we won t visit these cases or this method You need to take a course in ODEs although you could also selflearn from the text which is very clear on these subjects Let s take a closer look at the special case of Bessel s Equation which is well known and often encountered in hydrology From the text p 189 55 134 New Mexico Tech Hyd 510 Hydrology Progmm Quantitative Methods in Hydrology One of the most important ODES in applied mathematics in Bessel s equation6 l rzy w39 X2 7 V2 0 Its diverse applications range from electric fields to heat conduction and vibrations see Sec 129 It often appears when a problem shows cylindrical symmetry just as Legendre s equation may appear in cases of rphcrica symmetry The parameter 1 in l is a given number We assume that v is real Bessel s equation can be solved by the Method of Frobenius see text p 189 Accordingly we substitute the series 3c 2 ym 2 unim mg i 0 1470 with undetermined coefficients and its derivatives into 1 This gives x as 2 m rm r 7 lamx m quot E m ramrm quot 7710 7710 x 3o E ITH Jll l r lnz 7 V2 2 dmxmerr O 7110 1110 Etc leading to p 202 56 General Solution of Bessel39s Equation A general solution 3136559 395 cquurionfbr all values Q7quot 1 and r gt 0 is 9 39 C139 CZYIM39 7 lymxzwt 22mnml n m I 11 do x E m0 Jnx is called the Bessel function of the first kind qfordcr n The series l l converges for all x as the ratio test shows In fact it converges very rapidly because of he factorials in the denominator 135 New Mexico Tech Hyd 5 Hydrology Program Quznlitative Methods in Hydrology a mo me cos m 7 LAD 7 sm v17 1 Ynx lim Yyx This function is called the Bessel function of the second kind of order 1 or Neumann s function7 of order 1 Figure 109 shows YOU and Y xl fur n A n wu uhluln ll39ulu l I I l IIIL Bessel rumnon nl urdrr I x 7 l m zm 12 Jar V We 77 Al l A 2mhnhz 7 F il VhlL39ll limb mmlur II n cmmc mg lU7L For n l w olmun ll Bessel function of order 1 x 2m1 4 I Imx I 13 I m E o I 139 quot11 m II 1 m0 MIMI IUHL lllllllll In I lll39 lFlg lll7l Blll llIL crux ul llIL w r llh l l ll quot Ilnl L muplclcl n INC HIM Tuhlc Al quotI App SI and Ill Imglu of Ihc quotWain IlEIII39CllnEs Illl illclt llslllg Ilcurh lg X In U I A A l 7 4 L FII quotIquot Bessel funCtIons 0f ths first klnd ll and l x H quot hm hH Y r J m In 39 x3 quot 17 quot 2 y 77 IE 22 quotul m I39 II 8 s l 11 i l Z 2 quII u quot 7quot HI whch 39 r i u 0 ll unll as in 4 II I Ill 2 I I I l l I I l l l I I m 2 m quot39 L 2 IN H Y 06 quot Y I I I39I u 39 39 1 P 110 34 lt01 39 r h m Bessel functions of the second kind YU and Y For a small table see App 5 136 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Numerical Solution of 2quotd Order Linear ODEs We re still looking for solutions of the general 2quotd order linear ODE y pxy qxy rx with pq and r depending on the independent variable Numerical solutions can handle almost all varieties of these functions Numerical solutions to secondorder Initial Value IV problems can be solved by a variety of means including Euler and RungeKutta methods as discussed in 213 of your text We won t discuss these applications here as we don t have many 2quotd order IV problems in hydrology However we have lots of 2quotd order Boundary Value Problems BVPs In BVPs x usually represents space For these situations we use nite di erence methods which employ Taylor Series approximations again just like Euler methods for lSt order ODEs Other methods like the nite element see Celia and Gray 1992 nite volume and boundary integral element methods are also used The nite element method is the most common of these other methods in hydrology You may also encounter the socalled shooting method discussed in Chap 9 of Gilat and Subrarnaniam s 2008 textbook which you can safely ignore this semester As most hydrological BVPs are solved with the nite difference method that is where we ll focus our attention For example the popular groundwater ow code MODFLOW uses nite differences Recall the forward and backward exact Taylor expansions for yx The forward expansion is sz m3 2quot y x 3quot y xi 1 yxnl yxn Ax y39xn where step size Axquot xquot 1 7 xquot It projects forward in x from xquot to xn The backward expansion is this is shifted by one step compared to the way we wrote it before Ax 2 Ax 3 5 1 yquotxn 5 1 y39 xni 2 yxH yxnAxH y39xn We can manipulate these to get models for the rst and second derivatives in a 2quotd order BVP Your textbook discusses this on pp 910911 but in the context of PDEs The idea is just a generalization of the one we present here It is easier to do this if the step size is uniform Axquot Ax h is a constant We and your book make that assumption 137 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology 0 0 A locatlon x1 x2 xn1 xquot xn1 xN1 xN node 1 2 711 n n1 N l N Notice that I ve also relabled the rst node point to be x1 while your text insists on labeling the rst point x0 The subzero is a pain when indexing vectors and matrices in matlab and other programming languages For IV problems as in Text 211 this doesn t cause any problem since these values or the gradient are known and the solution propagates in one direction But our 2quot order problems are almost all BV problems and as we ll soon see we have to solve simultaneous solutions using matrix algebra Starting the index with the number 1 is very helpful as the rst row of the column vectors and the matrices that we derive will then correspond to the rst node point located at x1 the n row corresponds to the nth node point and so on Aside there is another way to handle this We could have still started node 1 at location x0 but labeled the node there as node 1 Then the last node on the right side of the domain would still be node N located at say xL In other words we can be sophisticated and de ne an indexing map between node numbers and node locations At the moment I prefer to take advantage of being able to equate location and node numbers Besides forward and backward differences how else can we approximate the rst derivative at a point xquot in a BVP Instead of using either 1 or 2 to build a model for the rst derivative we can use an average of both That is we subtract 2 from 1 and solve for the rst derivative This eventually leads to the socalled center dz39 rence numerical approximation of y 09 To get the second derivative expansion of y xn we can add 1 and 2 We ll do this below rst by looking at exact Taylor expansions and then their nite difference approximations From these results you can investigate the truncation error of the numerical approximations to each expansion Exact Taylor Series Expansions for 1St and 2m1 Derivatives and the ODE Rewriting l and 2 for constant step size we have yxyxAxy39xA yquotxA 3f y 39x 3 yxyx Axy39xA wow A y 39xi 4 where 3 is the forward expansion and 4 is the backward expansion and we ve assumed a constant step size Each of these can be solved for the derivative y xn as we previously did on pages 9799 of these notes The forward 6 and backward 15 equation numbers from pp 97 99 expansions for the derivatives were Ax 3 y 39xn 5 yxn1yx L Mar Ax Ax NZ H y 95 138 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology yx yx 1 1 M M3 v x n H x quotl x i 6 yn M M 2yn 3yn In the second equation we ve shifted the backward difference by one step to get derivatives at xquot We can also get a central expansion for y xn Subtracting 4 from 3 cancels out the even derivative terms and yields Ax3 Axs m 3 y xn 5 yxn1yxn12My x y xn 7 which we then inverse to solve for an exact value of y xn in terms of these two Taylor Series or m5 5 y39 xn 8 W yx yx g 3 m x 2Ax Ax 3 y n In 5 6 and 8 we have three exact expansions for the rst derivative at xquot If we truncate these expressions we get numerical approximations to the rst derivatives Note the relative truncation error that will result by comparing 8 to 5 and 6 The truncated central expansion is higher order than truncated forward or backward expansions For 2quotd order ODEs we need a similar expansion for the 2quotd derivative We get that by adding equations 3 and 4 which cancels out the odd derivative terms and yields sz Ax4 yxn1yxn12yxn2 2 yquotxn2 4 y quotxn 9 We solve for the second derivative to get H yxn12yxnyxn1 1 Ax y xn M m2 4 y xn 9 We now have expansions for both rst and second derivatives and can substitute them into the ODE We re looking at 2quotd order ODEs of the form y Pxy39Qxy V06 10 If we examine this equation at the point xquot in our domain we rewrite 10 as y xnPxny39xnWan9V06 11 139 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Then when we substitute our Taylor Series expansions we still have an exact expression The only decision at this point is which of the three expansions for the rst derivative to substitute In practice all three are used depending on the problem but most often it is either the backward 6 or central 8 expansion that is used The choice has to do with numerical truncation and stability issues that we encounter when truncating the expansions and solving the resulting nite difference equations numerically In commercial codes it is often a user decision so you need to be aware of this option Let s use the central expansion Then 11 becomes yxn1 2yx yxH Axl 6106quot yxn V06quot M M 12 mn23from We could also gather like terms together and rewrite this as yxn1 yxn71 x P n 2m 1 pm 2 1 pm 2 m2 2m yxn71mz Qxn yxnsz 2Ax yxn1 rxn ml 4 l l l l m2 l l l I y xn23y xn Note that besides the forcing term r the right hand side contains only higher order terms in the expansions These are the terms we ll neglect in our nite difference approximation and they define the local truncation error As you can see using centered expansions for both the rst and second derivatives results in higher order terms of order sz If the step size were to vary in x then the truncation error is not as good It s of order Ax 13 Finite Difference Approximations We can truncate all of our expansions and write nite difference approximations indicating the order of the approximation truncation error Then the exact value yxn is replaced by its numerical approximation yquot and the derivatives are replaced by their truncated expansions We did this once before for the forward and backward difference approximations to the rst derivative when discussing the Euler method In any case here are the results for the present case taken from truncating 5 6 8 and replacing the exact values of y with their numerical equivalents For the rst derivative approximations v N yn1 yn Forward dlfference y xn T 0Ax 14a Backward difference y xn E y yHOAx 14b 140 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Centered difference y39 xn E OAx2 140 For the second derivative approximaton 2 yquotxzW0M 15 For the 2quotd order ODE y px y qx y rx using the centered difference 140 for the y term we get from truncating 12 and 13 respectively 2 y my y pxgtyquot 2A mm rx 0W 16 01 1 pxn 2 1 pm 2 m2 2m yH Ax2 Qxnyn Ax2 2m yn1 rxnOAx 17 Some people prefer to write the difference version of the ODE as in 16 and others as in 17 There is also a difference in programming styles that probably leads to these biases as we ll see later when writing pseudocode Also many people prefer to multiply l6 and 17 by Ax It leads to a cleaner product We won t do that here There are N equations 16 or 17 one for each node When solving these equations we drop the truncation error terms terms of OAx which are only used theoretically to give us some idea of the error Also observe that each equation is an algebraic equation By invoking the nite difference approximations we ve reduced a single ODE to a system of N algebraic equations which are presumably easier to solve Since each node in the interior of the domain has two neighbors one to the right and one to the left each algebraic equation must be solved together with its neighbors This chains all of the nodes together and the N equations must be solved simultaneously It is important to grasp that this connection is omnidirectional unlike the Euler method for IV problems which propagates in one direction Because we have simultaneous algebraic equations we can use matrix algebra to solve them The terms in brackets on the left side of 17 become the entries along the diagonal middle term next lower diagonal left term and next upper diagonal right term of a tridz39agonal coef cient matrix A see below The forcing or load rxn on the right goes into a column load vector b The unknowns yquot go into a column vector y of unknowns New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology The N simultaneous equations 17 one equation for each node 711 2 Nl N are then written as Ay b 18 and solved for y and thus the nodal y s using matrix equation solvers like Gaussian Elimination LU factorization or some other direct matrix solution technique or by iteration we ll talk about this later Since there are N node points the load vector b is NXl b V1V2 rn1rnrn1 m m 19 where r n rxn The column vector of unknowns often called the state vector is also NX y Mal2 yn1ynyn1 yN1yNT 20 The square coef cient matrix A is NXN If 61 J is the i3911 row and j3911 column entry of the matrix A then the three nonzero diagonals along a row are arm1 anquot qxn 21b am 923 21c while all other enties aid are zero For an example 5node domain N5 the general version of matrix A is 34 0 22 an al2 0 0 y1 b1 all 6122 6123 0 0 y2 b2 0 6132 6133 6134 0 y3 7 b3 23 0 0 6143 6144 6145 y4 b4 0 0 0 6154 6155 y5 b5 142 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology In practice we alter nonzero values of 61 J and bj in the rst row or two rows and the last row or two rows to account for the boundary conditions one at x1 and the other at xN The nature of the alteration depends on the type of boundary at each end We ll see more about boundary conditions below Notice that if the ODE 10 has a rst derivative term ie p i 0 then the matrix A is not symmetric aid a for the rst two offdiagonals jz39 l If p0 then the matrix is symmetric aid a Equation solvers can and should take advantage of this symmetry less computer storage and computational effort Equation solvers can also take advantage of sparseness For example the tridiagonal form of A e g see 23 leads to an especially ef cient and oftused Gaussian Elimination equation solver called the Thomas Algorithm If we approximate the rst derivative with the backward difference 14b then 17 would become instead LC 2qxy rx0Ax2 24 This alters the entries in matrix A if p is nonzero and while different A remains nonsymmetric for this case arm1 anquot LC 22qxn 25b am 25c The nonzero p rstderivative term is used to represent advection in heat and mass transport problems while the second derivative is used to represent diffusion The zero order q term is used to represent mass or heat transfer to another media e g from water to sediments or 1st order decay Notice that in this application it is the presence of advection that leads to non symmetry of matrix A Example 1 For example consider the steadystate solute transport model with advection diffusion and decay aC a where Cx is concentration v is advective velocity of a uid D is a diffusion or dispersion coefficient and a is a rstorder decay constant I ve used partial derivatives for a reason that will be clear below Rewritten in standard form we have 143 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology C0 b with yC xx p vD odD and r 0 in the generic form y39 px y qxy rx If we ignore diffusion then the problem reduces to rst order While we could write a code to solve second order equation b numerically and then run it with D0 this would not be a good idea The 2quotd order ODE nite difference method derived here embraces smoothness but 151 order ODE advectiononly problems preserve sharp fronts The nite difference method in 17 or 24 smoothes those fronts whether diffusion is present or not Use a 1st order method instead like the Euler method What about problems with diffusion but that are still dominated by advection One then lives with the resulting nonrealistic numerical smoothing or adopts a better scheme This is one reason for choosing 24 over 17 If the ow is from left to right using the backward difference on the advection term helps to address this problem It s called upwinding or upstream weighting in more advanced methods If the ow were from righttoleft one would use a forward difference instead It s determined by how the grid is oriented to the ow The term upwinding implies one goes toward the direction from which the wind is blowing Another application is vadose zone ow It s gravity that is causing the nonzero p term while capillarity causes diffusion In this case upwind is upward While I ve raised this issue here it is mainly of importance in solving PDEs such as a with C x t and an additional solute storage accumulation term aCat on the left hand side see eqn b on p 120 of these notes Then one chooses between both spatial and temporal forward backward or centered difference schemes to manage the problem with tradeoffs with each combination If we ignore advection or there is none then p0 and the nite difference method in 17 or 24 works wonderfully This is the case with heat conduction solute diffusion and groundwater ow explaining the common use of the nite difference method to solve these problems Boundary Conditions On page 122 of these notes we talked about boundary conditions for 2quotd order BVPs Over a nite interval a 5 x 5 b there are two BCs one at each end of the interval In general the boundary conditions can be written as k1 ya kzy39a K 11 ND lzy b K where the k s and 1 s are parameters weighing the relative contribution of y and y at each boundary and the K s are values Recall the rst type or Dirichlet boundary For example if you know the value of y at a then k1 1 k2 0 and K is the value ya Then there is the second type or Neumann boundary Say if you know the gradient of y at b then 11 0 l 1 and Kb is the value y b Both k1 and k are nonzero at a third type or Cauchy boundary We can implement all three boundary condition types in the nite difference method We ll focus only on lSt and 2 type boundary conditions In practice for big problems it is preferred to apply 2n type BCs before lst type BCs However it is easier to learn by applying lSt type rst 26ab nd 144 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology KEY POINT Finally a point not made previously that applies to all solutions of 2quotd order ODEs you need at least one I or 3 type BC If you make both BVP BCs of 2quotd type you won t get a unique solution try it with a problem you can solve by the method of constant coefficients In numerical methods the computer will pick one for you but it will be arbitrary This can be embarrassing Dirichlet BCs There are several ways to handle lSt type boundary conditions The simplest and most elegant method but not the cheapest method for complex spatial domains and PDEs is to simply partition matrix A so that the nodes with known heads are no longer in the vector of unknowns and matrix A is reduced in size by the number of nodes with known values of y We have to explain this by example Suppose that at the rst node located at x1 the value of y is known Then we know exactly that yl yx1 and we don t have to compute it We can ignore the rst row of our 5node example in 23 However yz depends on yl and we have to respect that The second row of the computation in the 5node matrix equation 23 is az1y1 022 a23y3 b2 27 However we know yl We can move it to the right side creating a new entry in the load term for this 2quotd node point a22y1 6123 b2 a21y1 28 We partition the matrix A in our example 5node matrix by removing the rst row and rst column yielding a new smaller but still square matrix A22 or all a12 0 0 all I a12 0 0 0 a1 an an 0 Z ra2 ZZZ 6393quot A A A 0 as am as 0 a32 am as 0 A11 All 29 0 0 a43 a44 a45 0 0 a43 a44 a45 21 22 0 0 0 a a 0 0 0 a a 55 where 4x4 subpartition a a a 0 A22 32 33 34 0 6143 a 6145 0 0 6154 6155 We make similar adjustments to b b1 b2T where b1b1T and b2bza21y1 b3 b4 b5T Then the 5node matrix problem becomes 145 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology an an 0 0 yz b2 a21y1 Q32 Q33 6134 0 ys b3 31 0 6143 6144 6145 Mt b4 0 0 6154 ass y 5 b5 In general A22 yz b2 where in our example y2 yz y37y4 y5T is the vector of unknowns having removed them from the vector y which contained both unknowns and the known BC values of y Note that this approach works exactly the same way if the other xN BC is also rst type We would partition again ending up with a 3x3 partitioned matrix which we would still call matrix A22 and a new 3x1 load vector still called b2 but now b2 bza21y1 b3 b4a45y5T In fact this way of handling lSt type BCs applies to the spatial portion of PDE numerical solutions whether 1D 2D or 3D in space to any number of lst type nodes and to many different numerical methods e g nite difference nite volume nite element boundary integral element methods Newmann BCs There are also several ways to handle 2quotd type boundary conditions and their application depends on the method of derivation of the nite difference approach Our nite difference approach is called the point centered nite dz krence method The typical way to handle 2quotd type boundary conditions in this approach is to use ghost nodes These are pseudo nodes that lie outside the domain along 2quotd type BCs and are used to control the gradient y at the boundary In effect we assume that the same ODE applies outside of the domain across the second type boundary for a distance of Ax For example consider a 2quotd type zerogradient boundary at the node located at xN A sketch of this boundary and the solution might look like zero gradient boundary boundary ghost node domain To maintain a zero gradient for uniform grid spacing the value of y at the ghost node at xN1 is set equal to the value of y at the rst interior node located at xN1 or yN1 yN1 If the 2quotd type 146 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology BC involves a known nonzero gradient y xN we then apply the central difference 140 to preserve it yN1erl2Axyix1 33 As you can see it is just a generalization of the zero gradient boundary We now need to add a ghost node to our original problem Let s proceed with our 5node example where node 1 is a 1st type node that we ve already addressed The resulting matrix equation is 31 Now add a 2quotd type BC at xN of the type shown in 33 We do that by conceptually adding padding an additional row and column to the matrix and vectors leading to 22 6123 0 0 0 y2 bZ auy1 6133 6134 0 0 y3 b3 0 6143 6144 6145 0 y4 55 ys 6166 y6 34 a a u WWW a where the subscript 6Nl refers to the ghost node and according to 33 y6 y4 2Axy xN The values of a575 61575 and 61575 come from 21 or 25 depending on the choice of first order difference Ignore the last line Instead consider the second to last line representing computations at node N 5 6154y4 f assys f assys b5 35 Substitute the boundary expression 33 for y5 6154y4 a55y5 a55y4 2 6155 My OCN b5 36a or 6154 615514 assys b5 2 6156 AXy OCN 36b The problem reduces back to that of four unknowns y2 yz y37y4 y5T with a further revised load vector and a revised matrix an an 0 0 yz b2 a21y1 Q32 ass 6134 0 ya b3 37 0 6143 a 6145 Mt b4 0 0 6154 ass ass ys b5 26156 Axy xN Note that if the matrix A was symmetric before it is not now Also you may wonder where to get the value of 6155 since it appears to describe conditions just outside the domain But if you inspect 210 or 250 you ll see that it depends only on the step size and for the pterm central difference on the value of p at node 5 which is on the boundary and is de ned 147 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Aside There is another major approach to finite differences involving describing blocks rst and locating node points at the center of blocks This block centered nite dz krence approach is use by MODFLOW and most oilreservoir simulators This method treats 2quotd type BCs differently They don t need ghost nodes In that approach zerogradient 2quotd type BCs are natural boundary conditions That is without telling the code anything else it assumes that the default boundary is zero gradient This occurs at the edge of a gridblock The nearest node to the boundary is half a block away If you want to have a nite gradient across such a natural boundary you apply a sourcesink term at that node Such natural BCs are part of the algorithm and not something you have to program in The nearest thing to a natural boundary in a point centered scheme like that discussed here is a 1st type BC Pseudo Code The basic structure a nite difference code for solving 2quotd order ODEs by the nite difference method is INPUT The step size number of steps BCl type and value BC2 and parameters p q and r Note 1 the length of the domain is o en used instead of the step size or number of steps length step size X number of steps Ax N 1 Input two of these and calculate the third Note 2 since the presence of p950 can affect storage and solution a ag indicating the absence of that term or simply a ag set by noting that it is zero is usually included Some codes simply require that p0 BUILD A Fill or build the coef cient matrix A BUILD b Fill or build the load vector b SOLVE Solve Ayb for the vector y ANCILLARY Perform backsubstitution to nd uxes at lst type BCs check mass balance and other ancillary calculations OUTPUT Output y and uxes POSTPROCESS Plotyn n 1 to N and as needed uxes This is the basic organization of most nite difference codes A principle difference between codes is in building A and b and in the solution Many codes contain a library of solution methods later As for building A and b the two main approaches are illustrated by the two ways of writing the ODE difference approximation in 16 and in 17 Each would have its own different pseudo code The code related to 16 looks at each term in the equation and builds A and b in sequential steps A crude pseudo code for building A and b is 148 New Mexico Tech Hydrology Program Matrix Size Hyd 5 Quantitative Methods in Hydrology First you might note the presence of lSt type BCs and aim the nal matrix size accordingly For a onedimensional ODE the matrix is of size NXN with two 3rd type BCs N lXNl with one lSt type BC and N2XN2 with two lst type BCs You would note this and partition the matrix from the start Loop over rows Then row by row or node by node you would build in the contributions from the rst term on the left side ofl6 For example on row 71 you would add lAx2 to any and 2Ax2 to anyquot the diagonal Next you would cycle through the rows again and add in the contributions from the second term on the left of 16 Using an as the example for row n you would add pxn2Ax to any and 0 to the diagonal Next would come the third term on the left of 16 You would add qxn to the diagonal Finally you would add the loads rxn to the nth row of load vector b At this point you would have built the basic versions of A and b but without accounting for the BCs Loop over BCs types and locations Next you would adjust the appropriate rst or last row of A and b to account for 2quotd or 3rd type BCs Finally again you would insert the adjustment to the appropriate rst or last row of b to account for 1st type BCs While the main pseudo code separates building A and b in this version those operations merge when it comes to BCs although they don t have to This approach to building A and b is typical of blockcentered codes It is very powerful when dealing with complicated multidimensional coupled processes problems with complex boundaries Point centered codes typically use a ll method instead more along the lines ofl7 These codes are typically restricted to simpler problems A pseudo code for lling A and b might be Matrix Size First note the presence of lst type BCs and aim the nal matrix size accordingly For a onedimensional ODE the matrix is of size NXN with two 3rd type BCs NlXN l with one lst type BC and N2XN2 with 149 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology two lSt type BCs You would note this and partition the matrix from the start The indices used below are for this reduced matrix For all interior nodes 100p over the rows lling the diagonals of A and the entry of b Fill any from 21a Fill a from 21b This is o en done by de ning three vectors each Fill arm from 210 representing the entries in one of the diagonals nnl Fill bquot rxn At x1 adjust b1 and ifneeded all for the BC At xN adjust bNand if needed aN17N for the BC Example 2 Consider essentially horizontal steady groundwater ow in a leaky con ned aquifer modeled by d d K K ODE T 39 38 A dx B 31 Q lt a BCl d 0 BC2 7 H 38bc dx x70 XiL where x is hydraulic head Tx is spatially distributed transmissivity K B is a spatially distributed leakage coef cient 1 x is the known spatially varying head in the overlying phreatic aquifer and Q W x is spatially distributed pumping sign positive for injection negative for pumping BCl represents a no ow boundary BC2 represents a fullypenetrating known head of value H a lake Rewrite the ODE in standard form7 d2 idl OQWy abc2 T dx dx TB TB In terms that we are using in this course where the general ODE is written as y px y39qx y rx our unknown and parameters are x a 40a px d lnT dx 40b qx K TB and 400 W K TB o Q wT qx o Q wT 40d 7 In general it is not a good idea in finite difference methods to break up the 2nd derivative into a first and second order terms for spatially variable coefficient T as l have done here It is better and common to develop a 2nd order difference operator that explicitly handles this variation Such issues occur in a wide variety of diffusionlike problems Why Recall that 2quot derivative approximation of finite differences in space smear solutions and breaking it up this way leads to additional error Also note the derivative of T which itself represents a problem if T is not a continuous function of x 150 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology For this problem which is basically diffusive we should use the central difference approximations of 21 Thus for interior node points an 41a anquot Li mom 41b 61mm 410 b roe 41d where p q and r are de ned in 40 We can then use a ghost node at x1 and partition the matrix at xN respectively to represent BCl and BC2 Example 3 For simplicity lets assume Tconstant K B constant and do constant We ll let the pumping vary with x Then p0 q constant and rconstant fcn of x and the general equation is oftype y qy rx This simpli es the entries for the interior nodes 712 3 N 2 N 1 a1 am 42a anquot 51 42b bquot q 0 Q w xT 420 At node 1 we need to account for the no ow zerogradient boundary condition 6112 A2 43a b1new xl q o Q wJClT q o Q wmT 43b The last equation simpli es since the gradient is zero at the boundary At node N I the last node with an unknown head in this model we need to revise the load vector H q 0 Q w xNT 44 b5 new 39 1 ml where H is the known head at xN New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology Example 4 Let s simplify further by assuming a six node model including the known head H That leaves 61 5 unknown heads Let s also assume that the pumping only takes place from node 3 located at x 2m Finally let s multiply through by sz to get rid of those pesky denominators The interior node n2 3 4 5 contributions to A are am ann1 1 45a anquot 2 qsz 45b The load vector at n2 4 5 is bn q o M2 450 while at node n 3 bquot q 0 m2 Q wx3 szT 45d At node 1 6112 2 46a b1 q in M 46b while at node 5 b5Hq 0Ax2 460 where H is the known head at x5 Written out the matrix equation becomes 2qAx2 2 0 0 0 y1 qAx2 0 1 44sz 1 0 0 y qAx2 0 0 1 45sz 1 0 y3 qszg O QAxZT 0 0 1 2qAx2 1 y qAx2 0 0 0 0 1 44sz y qszg O H 47 where yquot n and Q Q w When p0 the resulting homogeneous equation y y 0 is sometimes called the Poisson equation This is a solution for that equation in 1D for the given BCs discretization and poin sink pumping Example 5 Suppose that there is no leakage so that q0 Then the ODE becomes a 1D version y 0 of the LaPlace equation sz0 in this case with a point sink at x 2m Then 47 simpli es to 152 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology 2 2 0 0 0 y1 0 1 2 1 0 0 yz 0 0 1 2 1 0 y3 QAx2T 48 0 0 1 2 1 y 0 0 0 0 1 2 ys H Note that the finite difference operator for the LaPlacian has the form 1 2 1 This is called the LaPlacian stencil and is found in basic nite difference solutions of diffusive mass transfer heat conduction and groundwater ow in a homogeneous aquifer Note how BCl the zero gradient boundary alters the stencil 153 New Mexico Tech Hyd 510 Hydrology Progmm Quantitative Methods in Hydrology Systems of ODEs Chapter 4 your textbook introduces systems of first order ODES In general these can be represented by the matrix expression y fty where y yl yz y3 y ynT is a column vector of unknows tis a scalar independent variable and the prime indicates differentiation wrt to t Typically for us the independent variable tis time This can also be written as shown below taken from p 134 42 of the text The first order systems in the last section were special cases of the more general system l i fi i an r fo Vls 39 gt Tn l r l n qu 39 i s We can write the system I as a vector equation by introducing the column vectors y D l 39 39LT and f f1 fur where T means transposition and saves us the space that would be needed for writing y and f as columns This gives u y fU y This system 1 includes almost all cases of practical interest For n 7 1 it becomes 1 ha 1 or simply I ft y well known to us from Chap 1 A solution of l on some interval a lt I lt b is a set of n differentiable functions Yl 11 quot39 In n on a S t S b that satisfy 1 on this interval In vector form yht hl hz hg hm kn An initial value problem for 1 consists ofl and n ICs y1t0 K1 y2t0 K2 y3t0 K3 ynt0 Kn where the K s are constants or ytoK KbKLKL Ki1KnT In all cases that you will see in hydrology systems of equations like that in l are IVPs after all it is a system of lst order ODEs Coupling Equation 1 represents n coupled equations which can be linear or nonlinear Reasons for coupling You are likely to run into two cases of coupled equations and their combination 154 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology In the first case the y s represent different but coupled processes For example consider uid ow in a porous media and coupled heat transport There are two coupled equations n2 yl would represent temperature and yz hydraulic head Add in a solute and you get a third equation since now n3 that represents solute concentration y3 The equations are coupled since both temperature and solute affect uid density and viscosity while uid ow advects heat and solute In the second case the y s represent discrete nodal values of a spatially distributed unknown such as concentration where the spatial domain that has been numerically approximated by nite differences nite elements or some other method where there are n such node points In this case a timespace PDE like time dependent advectiondiffusion has been reduced to a set of n time dependent ODEs one for each node and space no longer exists That is y remains continuous in time t but discretized in space Thus y23t would represent the value of y eg concentration at spatial node 23 at time t If 1 is solved numerically say via the backward Euler method the system of simultaneous ODEs becomes a system of simultaneous coupled algebraic equations that march forward in time Simultaneous or sequential coupling In practice coupling can be sequential or simultaneous Another set of equivalent terms is explicit or implicit coupling In simultaneous coupling equations 1 are solved together simultaneously or at least approximately so If 1 is linear this is done using linear matrix algebra with a matrix equation solver like Gaussian Elimination It is implicitly realized that each y depends on the others In sequential coupling one equation one line of 1 is solved and fed into the next and so on The assumption is that y depends on y1 but not the other way around This works perfectly well for systems which are truly sequentially coupled but is a very crude way and often inaccurate or even unstable way of handling problems that are actually simultaneously coupled Nevertheless it is sometimes employed in hydrology Consider as an example the common assumption that ow does not depend on temperature ie ignore the dependence of uid viscosity and density on temperature but heat transport and temperature depend on ow Solve the ow problem then use that solution to advect heat Mixed sequentialsimultaneous solutions are sometimes used For example a common approach to solving coupled multiphase ow oil water and gas is to use an IMPES scheme meaning implicit IMP pressure explicit ES saturation If it fails to work then most codes give you the choice to use a full simultaneous solution approach In any event your book focuses on true simultaneous coupling as will I Linear Systems text s 42 43 44 y 46 39 First carefully read Example 1 in 4l on p 130 of the text before continuing It Wmmm 9 an an describes a linear problem with two coupled mg 77 serum mm mTanksTlowerand r 7 equations describing solute mixing between two tanks This has applications in hydrology for example mixing of a tracer in two tanks prior to injection in an aquifer or stream More commonly this is analogous to mixing between two hydrologic reservoirs such as wind driven circulation mixing of two portions of a lake 155 Nstexlcu Tech HdelEI Hydrulugmegam QuanmauveMmhuds m Hydrulugy Mixing Problem Invulvlng Two Tlnks ml pmblcm umxlung mglv Link I mauled by mm um um IwH um 4 mm m cnnewpummg Eumplc m x w 1 beeJust m pnnmple of modeling um be m Ame W um um Th quotmun wquot Iv u mm m um nmmmcr an a T k1quot nn T mum gum Hum v 14anka Icreaclh In 7 m we pun thruM 150 w r r l D Alml T WIquot umum a lea mu In much lemllzcr as mm IH h Inn m Ti 39 Salun un Step I 39 E A fur A u y L 39 l q m xm mmm numnu Swmxhul an un T lmm n 77 m m m l 7 ln nwmm 7 Om uwmul 7 Wk 7 1 ATunk r u hmnwImn nummvmm 7 7W W x 1a Helm me nwhcmicul mch u our mmm pmlgkm n m mmn m hnlmnlur mm 34 7mm A 01m hqu I 39 uozy 7 002 mm 121 m gt l wumm m We need quotmumin umunm m unnm mm m um ldznncnlly mm chm m am 1 me rnr eigenvalue 00 7 A ltan m an m 7 A 7 nu 7 m2 7 mqu 7 m 4 mm7 u rm 7 n 02 7 Mm an W g mu A2 70 04 Eigemuun um ubminuu from 114 m 521240 mm a u and A 7 7004 quurprcwnx A mi gives lwc nccd uuly mu nm equulmn m x m 7mm mm 7 n and van 4 com 77 cm2 7 mpmwely Hm r and 7 7 1 mpmwely md w can Ink x a 7 1 arm 7x 7 v om A 7mm y 1 v 0Wquot wow 62 0mm mm NewMexieo Tech Hyd5 Hydrology Program Quantitative Methods in Hydrology From n ma mt uipcnmnmn pnmpic Mmn eummuek m hnhl for Nam u luvnmgrnrnm linuur H w mm quotmm A oluuun 4 l i y l mg A ma a l 27 l 7 a 7 m when q and 3 im rhumn commie mer we shall call us a general sol Iinlx ul39inhin condillulbi Thu minul umdmum m Him 0 mt lcmlucr in mm warm with 50 l l up 1 ls I39nnu lIi md 1 HI I 0 WI nb n n willA 2lllliif lrill Incuuumucuu Ilm pl 2 7 n r 7 2 my Th1 sulunun ml 7 7 we 7 7 75 11m plvm the answer 7 73 In rumpum nl h 7 75 75quot39 quot mini 7 Imm runel v2 7 75 NFMH39 unk 739 upper t39urvl39i mm 7 1mm um cwmwmi amtu m m and mu Exponential titmm m r in me umumou Imm 75 lb u i a i V r 7 nu cun cslnuk quotsymmetricquot wlum ml min clungc ll rl mlllnll cnnmmx39d mo in m fcmlncr and 7 ummmcd 50 lbquot Slfp 4 Answer 139 mnmms lull me icmlncr umhum m r n u mmums 3 11 he mxul mnmInL Ihul Ix 5n n I39Im 7i 7 75 quotM quotM g 7 lln Emmi 7 275 Mean mi rluul himld uruuim m u lcusl ll oul Imlr 1m huuv 7 15in Let s explore the mathematics behind this solution The system of equations in 1 on p 154 ofthese notes is y fty It represents any system of 1539 order ODEs If however the system of equations is linear it can be simpli ed to y Ay g as New Mexico Tech Hyd 510 Hydrology Progmm Quantitative Methods in Hydrology described below and as taken from the text p 137 and 46 In this last section 01 Chap 4 we discuss methods for solving mnhomogeneous linear systems 01 ODEs l f Ay g see See12 load where thMeelor gm is not identically zero We assume gm and the entries 0139 the n X 11 matrix AU to be continuous on some interval J of the I axis From a general solution y t of the homogeneous system y39 Ay on J and a particular solution y quot 01 l on J ie a solution of 1 containing no arbitrary constants we get a solution 01 l In 2 y yth y y is called a general solution at l on J because it includes every solution 01 t l on J The linear matrix equation 1 can be rewritten as a set of simultaneous linear algebraic equations For example if there are two unknowns y1 and yZQ then these equations are I I 1 171111 d121392 g1 1 51 1 2 2 2 I for example 1 1 2 21 1 2212 g2 quot2 13h 1 2 1 The full linear system of n firstorder ODEs in 11 unknown functions hm mm is of the form I u 7quot l The entries ay s represent 1 11H lzdl 1 1 g1 links between two dependent V7 7 variables y and y Whennis 2 quot211 quot2212 quotan g2 large A is usually very spam as eachms connected or coupled to relatively few of y the other States y Yn anlyl 122 2 39 39 39 11ml gr ASIDE Recall the 2nd order ODE BVP nite difference numerical solution where we wrote Ayb where b was the load vector The sign on either the coefficient matrix or the load vector must be changed to be consistent with the signs in Chapter 4 of the text Or just let g b Matrix Vector Notation We ve previously introduced and used matrixvector notation In any event our text reviews this in 40 and gets into more detail in Chapter 7 However our previous discussion in class was limited to matrix algebra We needed to de ne y in 1 p 125 and 126 of the text to get the result above Page 127 of our text extends the discussion to include 158 NeWMexico Tech Hyd510 Hydrology Program Quantitative Methods in Hydrology Differentiation The derivative of a matrix or vector with variable entries or components is obtained by differentiating each entry or component Thus if no 72 yl tr 72r2t ytr then y39r 2i sin i t cos I Using matrix multiplication and differentiation we can now write t l for g0 and n2 as Ni 11 12 1 5 2 1 7 y Ay em y v 21 22 2 13 1 2 NIH Superposition Principle We set g0 in the linear equation 1 on p 126 to get the homogenous equation y Ay Superposition Principle or Linearity Principle Ifym and ya are solutions of the homogeneous linear system 4 on some inicn39ul 50 ii any linear roziibiimiimi y Cly l39 c zy zl Homogeneous Solution text 1 137139 When the text refers to equation 4 it is referring to the homogeneous equation 4 y Ay From 1 138 the text By a basis or a fundamental system ol solulions ol the homogeneous system 4 on some interval J we mean a linearly independent set of II solutions ym 39 y of 4 on that intervall We write J because we need I to denote the unit matrix We call a eorrespomling linear combination 5 y ply 39 Cnym c1 39 1 arbitrary a general solution 014 on J 11 can be shown that it the ujkti in 4 are continuous on J then 4 has a basis of solutions on J hence a general solution which includes every solution ol 4 on J We can write it solutions ya ya ya yquot391 yquot of 4 on some intervalJ as columns of an moi matrix 6 Y y y y y y quot 159 New Mexico Tech Hyd 5 Hydrology Progmm Quantitative Methods in Hydrology The determinant of Y is called the Wronskian of ya y written rt xi rim 1 2 tn 2 2 39 2 7 Myth I you r W t1 2 n l n l n 39 77 The columns are these solutions each in terms of components These solutions form a basis on fit and only if W is not zero at any 11 in this interval Weither is identically zero or is nowhere zero in J This is similar to Secs 26 and 3i If the solutions ym 39 ym in 5 form a basis a fundamental system then 6 is often called a fundamental matrix Introducing a column vector c cl 72 39 MT we can now write 5 simply as 3 y Yc Method ofConstant Coef cients and Eigenvalue Problems 43 Reconsider the linear homogeneous matrix ODE equation 1 y Ay where tis the independent variable By de nition it has constant coef cients iff the a entries do not depend on t they already don t depend on y as the system is linear A single ODE of the form y ay has a solution of the form yC 9 This suggests a form for the solution of the matrix equation 2 y x e where x is a vector of constants and A is a scalar decay or growth coefficient of some kind Substitute 2 into 1 to get y lxel AyAxel Divide by c you can t divide by X it s a vector to get 3 Ax Ax If we solve this problem for A and x we can substitute into 2 and have the solution for l 160 New Mexico Tech Hyd510 Hydrology Program Quantitative Methods in Hydrology For every 1 there is a trivial solution x0 a vector of n zeros But there are 1 values for which the solution is nontrivial x 0 These 1 values are called eigenvalues of A and corresponding to each is a vector x called the eigenvector of the eigenvalue 1 Equation 3 is an eigenvalue problem see text p 129 There are n linearly independent eigenvectors x 312 xquot and corresponding eigenvalues 11 12 quot General Solution If the CUI ISIUHT nmtrix A in the system 1 has It linearly independent set of n eigenvectors then the Corresymming solutions ym y in 4 form a basis of solutions 0fl and the correspomling general solution is 1 t A f 5 y x 1 anx me quot ASIDE when solving timespace PDEs some numerical schemes discreteize in space n nodes points and then solve in continuous time by assuming a solution of this form Some elds call this the matrix exponential approach It is sometimes used in hydrology especially in socalled data assimilation schemes We can rewrite the eigenvalue problem 3 as A AIx 0 For x to be nontrivial the determinant of the coef cient A 11 must be zero see Cramer s Theorem of linear algebra inext page of these notes and its proof pp 312314 of the text If n 2 then we can solve directly using linear algebra For n gt 2 things quickly become messy as a studentI once required to solved an eigenvalue problem by hand with n5 it was tedious There are very efficient numerical solvers for eigenvalue problems see text s 206209 including those in Matlab For n 2 the expression A Mx 0 becomes see text p 129 11 T A 12 a 7 0 11 T1x1 T 12x2 0 7 or a 21 22 T 1 x2 0 21x1 T 22 T lx2 0 In this case the determinate of A M is 7 2 Det A 201 A M 1 a 12 which we set to zero or 12 T 11 T 22 1 T 11 22 T 12 21 0 b New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology This quadratic equation is called the characteristic equation of matrix A It has two roots for a solution which provide the two eigenvalues 11 and 12 for this two unknown x1 x2 problem Then for each root substitute it into simultaneous equations a and solve for the corresponding eigenvector perhaps by Gaussian Elimination yielding respectively x and 212 Your textbook has several examples Cramer s Rule pp 312314 of text Cramer s Theorem Solution of Linear Systems by Determinants a If a linear system of n equations in the same number of unknowns x1 n auxl auxz ulnx b1 21X1 22X2 agnxn 172 6 anlxl awxz t amxn b7L has a nonzero eoeyf cient determinant D detA the system has precisely one solution This solution is given by the formulas 7 rl g 392 X D Cramer s rule D D TL D where DR is the determinant obtained from D by replacing in D the kth column by the column with the entries 71 bn 1 Hence if the system 6 is homogeneous and D i 0 it has only the trivial solution x1 0 x2 0 nontrivial solutions x D 0 the homogeneous system also has Eigenvalue Problem Fiml the eigenvalues and ulgcllwcmn m the Innlrh 7H 40 11m A n l v l2 Solution The characteristic equation is llw quudmllc culminin 39 l A J w alli A All 7 l39 M t LI 4 n 7 I r L A ll link It whilium A z 3 and A n x 39Ilmeim1IwmpcmunuuI Fineanurx re uhmiued lmm 14 i Fol A Al 7 7 M h num ll l l t 4 v Inn 1411 s n Im 39HZ rllllxJill A mluutm m39 the ilrs cquulmn m l i 1 l2 7 llm in ulMlL39s he secuml euumxun WW i H un m clgcmctlvr Hi A mimpundmg m A 7 n i 1 I 1 M x Sumlull 39 7 I n b i m igumuclui nt39 4 wrmwmliu in A1 ll llvl1ylllvud rum l H n unh A A wmy lllh I 162 New Mexico Tech Hyd510 Hydrology Program Quantitative Methods in Hydrology Method ofConstant Coef cients Phase Plane Critical Points Stability Attractors 43 and 44 The later information in these sections on linear equations and in the next section 45 on non linear equations talks about solutions to coupled problems and introduces phase portraits It is important but it is beyond the scope of our effort this semester leu ni rnn minitun minim rnrm irnjmtnrtm m un phin plnlw nin n nn 39ul 4 rt innncxnninc t M and un uptiull h It TNIQEKOHFM Di the sy em Ill Vll quot Tlathom s 0V Ute yulzm llll Saddlu point lCautev This is the route to looking at chaos theory and limits to prediction how far into the future is a prediction reliable limited by coupling and nonlinearity in weather climate and hydrologic systems You start outwith simple coupled nonlinear equations and they reveal 50 of the issues encountered in large models If you are interested you can follow up on your own or take a math class I can also suggest additional reading for the future or we can also talk about it off line Method ofUndetermined Coef cients particular solutions 46 Recall from notes p 156 or text p 159 In llux LN cclinn nl Imp 4 we indn Illcllimk I whiny nunlmmngcncmh innn 39IL39III n1 nus IUc Sec 42I Item lIiL ct lm gm n Ilnl Identicth mm Wu imnnic gilt and the unitnu til the n nnnm m n lic continuum on some mlmull J in the l ann 1 pcimul wluliun y ll n the tnnnnnammn mmn y39 t on J nnl tcnlnr solullrm quotWII ul It nn I Lin t wlullun in II cmmuninp m nrlntrnry mnqnmq we gx l n wtuunn ul Ill 39239 in yr 3 ix uullctl l gumml sululion til I l i It I humth il lHL39lllI owl Union at him J The textbook continues on p 160 referring to the method of undetermined coe icients as an approach to the particular solution As for a single ODE this method uitable if the entries of A are constants and the components of g are constants positive integer powers of I exponential functions or cosines and Sines In such a case a particular solution yfp is assumed in a form similar to g for instance y p 11 VI Wt2 if g has components quadratic in t with u V w to be determined by substitution into 1 This is similar to Sec 27 except for the Modification Rule It suffices to show this by an example New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology See text p 160 for example We could also use the method of variation of parameters but we skipped that for a single unknown and will skip it again here for multiple unknowns Numerical Methods in independent variable t not in text Suppose we want to solve the nonhomogeneous linear model 1 y Ayg where tis the scalar independent variable A is a coefficient matrix not depending on y and y is the column vector of unknowns using a numerical method in t We simply scale up the lst order methods introduced before Euler methods For example we can use any QEuler method although it is usually a forward backward or midpoint Euler scheme that is employed Let s write the QEuler method from 37 on p 103 of the notes changing notation to that of l 2 I Ark6Akyksumac 6mmArk1 6g6gk where I ve replaced scale p on p 3 with matrix 7A the number 1 with the identity matrix I the load r with the load vector g the step size Ax with the time step At the step index n with the time step index k Notice that we ve maintained a variable step size Atk tk1 7 tk and that we ve allowed A and g to also vary in time It is easy to simplify 2 for either a constant time step or to a constantintime A and g Please note that all of the Taylor Series work and other derivation we did for the single ODE applies here especially truncation and stability analyses Aside Recalling that work you should also be able to relate 2 to the matrix exponential approach of 5 in 43 and p 159 ofthese notes and see how 2 above represents a truncated series expansions of the matrix exponentials and that the eigenvalues are essentially like the reciprocal of a time constant 7 actually 71 time constants where n is the number of eigenvalues Now you will see how the choice of weight 9 in uences coding as well as results If 90 as it is in an explicitforward in time FIT Euler method we can easily invert the coefficient matrix on the LHS as I391 I The result is 3 FIT scheme ym EIAtk AkykAtk gk As you can see we don t need a matrix equation solver The FIT or explicit method only involves inexpensive matrixvector multiplication Notice too that having A and g vary in time doesn t change the matrix math In fact A and g can vary with y as well as t and the problem can 164 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology be nonlinear The explicit method can handle it However the method is only conditionally stable as we discussed on p 106 of these notes The time step At is limited by the smallest of the n time constants of the system of equations that is the largest of the n eigenvectors 1 Let s consider the backwardintime BIT scheme with 91 Then 1 becomes 4 BIT SCheme 1 Ark Ak1yk1 E Yk Ark gk1 We have to invert the matrix I Ark Am a much more expensive operation especially as 71 gets large Actually we don t invert the matrix Instead we solve 4 by an equation solver like Gaussian Elimination depending on the algorithm this has from n2 to 713 operations Notice in 4 that if A and time step are constant in time your equation solver only need to decompose the coefficient 1 At A once say with LU decomposition You can then solve 4 at each time step by back substitution very inexpensively This is typical of tradeoffs in numerical solutions There is also the popular midpoint p 104 of the notes or CrankNicholson scheme 905 particularly popular when the y s represent spatially distributed node points of a spacetime PDE I ll leave that one as an exercise for you to write out You ve already learned how to build A for a nite difference approximation to a 2quotd order equation ODE in space Thus you should now be able to set up 3 or 4 for a nite difference method in space and respectively a FIT or BIT method in time From separate presentations on matrix methods you also know how to solve the matrix problems in 3 and 4 Without having written a PDE you are already set up to solve one Example Let s revisit Example 1 from p 140 of the text repeated above on p 156 of these notes Recall that it involves two unkowns for solute concentration in each tank represented by the homogeneous equation y Ay where y is 2x1 vector of unknowns and A is 2x2 coefficient matrix 002 002 A 002 002 Since the equation is homogeneous g 0 Thus the BIT Euler method 4 becomes 1 Ark AYk1 E Yk where the subscript on A has been dropped since A remains constant We can solve this example by the BIT FDM by Gaussian Elimination every time step Try it and compare to the exact answer It may help to write it out fully 1 0 002 002 y1 y1 Ark 0 1 002 002 yz M yz k or if you prefer 165 New Mexico Tech Hyd 5 Hydrology Program Quantitative Methods in Hydrology yum Ark 03902y1k1 03902y2k1 y1k yum Atk03902y1k1 03902y2k1 y2k Gathering terms 1 002Atky1k1 002Atky2k1 y1k 002Atky1k1 1 002Atky2k1 y2k or back to matrix form 1 A5002 A5002 yw1 ylyk A5002 1 A5002 ym1 y where it is this latest matrix that you assemble and then use in the Gauss Elimination solution for each time step Other methods in time You can also use other numerical methods in time While Euler methods are the most popular some codes use multistep methods For various reasons RungeKutta methods are not used to solve 1 for most typical hydrologic models unless y represents the two or three dimensional spatial coordinates of particle being tracked Euler Methods are also used In that case one is solving the more general problem y fty where f is the velocity vector at time t RungeKutta Fehlberg p 893 of text is popular for this purpose 166
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'