Systems of ODEs

1 New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 154 Systems of ODEs Chapter 4 your textbook introduces Systems of first order ODES. In general, these can be represented by the matrix expression y =f(t,y), where y = {y1 , y2, y3, .., yn-1, yn}T is a column vector of unknows, t is a scalar independent variable, and the prime indicates differentiation wrt to t. Typically for us the independent variable t is time. This can also be written as shown below, taken from p. 134 ( ) of the text: on a t b that satisfy (1) on this interval. In vector form y=h(t) = {h1 , h2, h3, .., hn-1, hn}T. An initial value problem for (1) consists of (1) and n ICs y1(t0)= K1, y2(t0)= K2, y3(t0)= K3, , .., yn(t0)= Kn, where the K s are constants, or y(t0)=K = {K1, K2, K 3, .., K n-1, K n}T . In all cases that you will see in hydrology, Systems of equations, like that in (1) are IVPs (after all, it is a system of 1st order ODEs).

2 Coupling Equation (1) represents n coupled equations, which can be linear or non-linear. Reasons for coupling. You are likely to run into two cases of coupled equations, and their combination. New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 155 In the first case the y s represent different, but coupled processes. For example, consider fluid flow in a porous media and coupled heat transport. There are two coupled equations (n=2). y1 would represent temperature and y2 hydraulic head. Add in a solute and you get a third equation (since now n=3) that represents solute concentration, y3. The equations are coupled since both temperature and solute affect fluid density and viscosity, while fluid flow 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 finite differences, finite elements or some other method, where there are n such node points.

3 In this case a time-space PDE (like time dependent advection-diffusion) 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 y23(t) would represent the value of y ( , 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.

4 In sequential coupling one equation (one line of (1)) is solved and fed into the next, and so on. The assumption is that y2 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 flow does not depend on temperature ( , ignore the dependence of fluid viscosity and density on temperature), but heat transport and temperature depend on flow. Solve the flow problem then use that solution to advect heat. Mixed sequential/simultaneous solutions are sometimes used. For example, a common approach to solving coupled multiphase flow (oil, water and gas) is to use an IMPES scheme, meaning implicit (IMP) pressure, explicit (ES) saturation.

5 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 , , , ) First, carefully read Example 1 in on p. 130 of the text, before continuing. It describes a linear problem with two coupled 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. New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 156 Or eqn. (a) on p. 162 of these Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 157 Let s explore the mathematics behind this solution.

6 The system of equations in (1) on p. 154 of these notes is y =f(t,y). It represents any system of 1st order ODEs. If, however, the system of equations is linear it can be simplified to y =Ay + g, as Questions: New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 158described below, and as taken from the text p. 137 and .. The linear matrix equation (1) can be rewritten as a set of simultaneous (linear) algebraic equations. For example, if there are two unknowns, y1(t) and y2(t), then these equations are [ASIDE: Recall the 2nd order ODE BVP finite difference numerical solution where we wrote, Ay=b, 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 matrix-vector notation.

7 In any event, our text reviews this in and gets into more detail in Chapter 7. However, our previous discussion in class was limited to matrix algebra. We needed to define 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: +g1+g2+g1+g2+gn+2 +1 The full The entries aij s represent links between two dependent variables, yi and When n is large, A is usually very sparse, as each yi is connected or coupled to relatively few of the other states yj (j i). load New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 159 Superposition Principle We set g=0 in the linear equation (1), on p. 126, to get the homogenous equation y =Ay. Homogeneous Solution (text, p. 137-139) When the text refers to equation (4) it is referring to the homogeneous equation (4) y =Ay From p. 138 the text: We can write n solutions y(1), y(2), y(2).

8 Y(n-1), y(n) of (4) on some interval J as columns of an nxn matrix, Y (6) Y = { y(1), y(2), y(2), .. y(n-1), y(n) } for g=0 and n=2, asDifferentiation. New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 160 Method of Constant-Coefficients and Eigenvalue Problems Reconsider the linear, homogeneous matrix ODE equation (1) y =Ay where t is the independent variable. By definition it has constant coefficients iff the ai,j 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 y=C eat. This suggests a form for the solution of the matrix equation. (2) y =x e t where x is a vector of constants, and is a scalar decay or growth coefficient of some kind. Substitute (2) into (1) to get y = x e t = Ay = A x e t Divide by e t (you can t divide by x, it s a vector) to get (3) Ax = x If we solve this problem for and x we can substitute into (2) and have the solution for (1).

9 W=|Y| New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 161 For every there is a trivial solution, x=0 (a vector of n zeros). But there are values for which the solution is non-trivial, x 0. These values are called eigenvalues of A, and corresponding to each is a vector x called the eigenvector of the eigenvalue . Equation (3) is an eigenvalue problem (see text, p. 129). There are n linearly independent eigenvectors x(1), x(2), .. x(n), and corresponding eigenvalues (1), (2), .. (n). ASIDE: when solving time-space PDEs, some numerical schemes discreteize in space (n nodes points) and then solve in continuous time by assuming a solution of this form. Some fields call this the matrix exponential approach. It is sometimes used in hydrology, especially in so-called data assimilation schemes. We can rewrite the eigenvalue problem (3) as (A- I)x = 0 For x to be non-trivial the determinant of the coefficient (A- I) must be zero (see Cramer s Theorem of linear algebra next page of these notes- and its proof, pp.)

10 312-314 of the text). If n = 2, then we can solve directly using linear algebra. For n > 2 things quickly become messy (as a student I once required to solved an eigenvalue problem by hand with n=5; it was tedious.). There are very efficient numerical solvers for eigenvalue problems (see text, s ), including those in Matlab. For n = 2 the expression (A- I)x= 0, becomes (see text p. 129) = 002122211211xxaaaa or 0)(0)(222121212111= +=+ xaxaxaxa (a) In this case the determinate of (A- I) is Det (A- I)=| A- I |= 22211211aaaa which we set to zero, or 2112221122112)(aaaaaa ++ =0 (b) New Mexico Tech Hyd 510 Hydrology Program Quantitative Methods in Hydrology 162 This quadratic equation is called the characteristic equation of matrix A. It has two roots for a solution, which provide the two eigenvalues, (1) and (2), for this two unknown (x1, x2) problem.