.. !split Many mathematical models involve :math:`m+1` unknown functions governed by a system of :math:`m+1` differential equations. In abstract form we may denote the unknowns by :math:`u^{(0)},\ldots, u^{(m)}` and write the governing equations as .. math:: \mathcal{L}_0(u^{(0)},\ldots,u^{(m)}) &= 0,\\ &\vdots\\ \mathcal{L}_{m}(u^{(0)},\ldots,u^{(m)}) &= 0, where :math:`\mathcal{L}_i` is some differential operator defining differential equation number :math:`i`. .. _fem:sys:vform: Variational forms ================= There are basically two ways of formulating a variational form for a system of differential equations. The first method treats each equation independently as a scalar equation, while the other method views the total system as a vector equation with a vector function as unknown. Sequence of scalar PDEs formulation ----------------------------------- Let us start with the approach that treats one equation at a time. We multiply equation number :math:`i` by some test function :math:`v^{(i)}\in V^{(i)}` and integrate over the domain: .. _Eq:fem:sys:vform:1by1a: .. math:: \tag{1} \int_\Omega \mathcal{L}^{(0)}(u^{(0)},\ldots,u^{(m)}) v^{(0)}{\, \mathrm{d}x} = 0, .. _Eq:_auto1: .. math:: \tag{2} \vdots .. _Eq:fem:sys:vform:1by1b: .. math:: \tag{3} \int_\Omega \mathcal{L}^{(m)}(u^{(0)},\ldots,u^{(m)}) v^{(m)}{\, \mathrm{d}x} = 0 {\thinspace .} Terms with second-order derivatives may be integrated by parts, with Neumann conditions inserted in boundary integrals. Let .. math:: V^{(i)} = \hbox{span}\{{\psi}_0^{(i)},\ldots,{\psi}_{N_i}^{(i)}\}, such that .. math:: u^{(i)} = B^{(i)}(\boldsymbol{x}) + \sum_{j=0}^{N_i} c_j^{(i)} {\psi}_j^{(i)}(\boldsymbol{x}), where :math:`B^{(i)}` is a boundary function to handle nonzero Dirichlet conditions. Observe that different unknowns live in different spaces with different basis functions and numbers of degrees of freedom. From the :math:`m` equations in the variational forms we can derive :math:`m` coupled systems of algebraic equations for the :math:`\Pi_{i=0}^{m} N_i` unknown coefficients :math:`c_j^{(i)}`, :math:`j=0,\ldots,N_i`, :math:`i=0,\ldots,m`. Vector PDE formulation ---------------------- The alternative method for deriving a variational form for a system of differential equations introduces a vector of unknown functions .. math:: \boldsymbol{u} = (u^{(0)},\ldots,u^{(m)}), a vector of test functions .. math:: \boldsymbol{v} = (u^{(0)},\ldots,u^{(m)}), with .. math:: \boldsymbol{u}, \boldsymbol{v} \in \boldsymbol{V} = V^{(0)}\times \cdots \times V^{(m)} {\thinspace .} With nonzero Dirichlet conditions, we have a vector :math:`\boldsymbol{B} = (B^{(0)},\ldots,B^{(m)})` with boundary functions and then it is :math:`\boldsymbol{u} - \boldsymbol{B}` that lies in :math:`\boldsymbol{V}`, not :math:`\boldsymbol{u}` itself. The governing system of differential equations is written .. math:: \boldsymbol{\mathcal{L}}(\boldsymbol{u} ) = 0, where .. math:: \boldsymbol{\mathcal{L}}(\boldsymbol{u} ) = (\mathcal{L}^{(0)}(\boldsymbol{u}),\ldots, \mathcal{L}^{(m)}(\boldsymbol{u})) {\thinspace .} The variational form is derived by taking the inner product of the vector of equations and the test function vector: .. _Eq:fem:sys:vform:inner: .. math:: \tag{4} \int_\Omega \boldsymbol{\mathcal{L}}(\boldsymbol{u} )\cdot\boldsymbol{v} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V}{\thinspace .} Observe that :ref:`(4) ` is one scalar equation. To derive systems of algebraic equations for the unknown coefficients in the expansions of the unknown functions, one chooses :math:`m` linearly independent :math:`\boldsymbol{v}` vectors to generate :math:`m` independent variational forms from :ref:`(4) `. The particular choice :math:`\boldsymbol{v} = (v^{(0)},0,\ldots,0)` recovers :ref:`(1) `, :math:`\boldsymbol{v} = (0,\ldots,0,v^{(m)}` recovers :ref:`(3) `, and :math:`\boldsymbol{v} = (0,\ldots,0,v^{(i)},0,\ldots,0)` recovers the variational form number :math:`i`, :math:`\int_\Omega \mathcal{L}^{(i)} v^{(i)}{\, \mathrm{d}x} =0`, in :ref:`(1) `-:ref:`(3) `. .. _fem:sys:uT:ex: A worked example ================ We now consider a specific system of two partial differential equations in two space dimensions: .. _Eq:fem:sys:wT:ex:weq: .. math:: \tag{5} \mu \nabla^2 w = -\beta, .. _Eq:fem:sys:wT:ex:Teq: .. math:: \tag{6} \kappa\nabla^2 T = - \mu ||\nabla w||^2 {\thinspace .} The unknown functions :math:`w(x,y)` and :math:`T(x,y)` are defined in a domain :math:`\Omega`, while :math:`\mu`, :math:`\beta`, and :math:`\kappa` are given constants. The norm in :ref:`(6) ` is the standard Euclidean norm: .. math:: ||\nabla w||^2 = \nabla w\cdot\nabla w = w_x^2 + w_y^2 {\thinspace .} The boundary conditions associated with :ref:`(5) `-:ref:`(6) ` are :math:`w=0` on :math:`\partial\Omega` and :math:`T=T_0` on :math:`\partial\Omega`. Each of the equations :ref:`(5) ` and :ref:`(6) ` needs one condition at each point on the boundary. The system :ref:`(5) `-:ref:`(6) ` arises from fluid flow in a straight pipe, with the :math:`z` axis in the direction of the pipe. The domain :math:`\Omega` is a cross section of the pipe, :math:`w` is the velocity in the :math:`z` direction, :math:`\mu` is the viscosity of the fluid, :math:`\beta` is the pressure gradient along the pipe, :math:`T` is the temperature, and :math:`\kappa` is the heat conduction coefficient of the fluid. The equation :ref:`(5) ` comes from the Navier-Stokes equations, and :ref:`(6) ` follows from the energy equation. The term :math:`- \mu ||\nabla w||^2` models heating of the fluid due to internal friction. Observe that the system :ref:`(5) `-:ref:`(6) ` has only a one-way coupling: :math:`T` depends on :math:`w`, but :math:`w` does not depend on :math:`T`, because we can solve :ref:`(5) ` with respect to :math:`w` and then :ref:`(6) ` with respect to :math:`T`. Some may argue that this is not a real system of PDEs, but just two scalar PDEs. Nevertheless, the one-way coupling is convenient when comparing different variational forms and different implementations. Identical function spaces for the unknowns ========================================== Let us first apply the same function space :math:`V` for :math:`w` and :math:`T` (or more precisely, :math:`w\in V` and :math:`T-T_0 \in V`). With .. math:: V = \hbox{span}\{{\psi}_0(x,y),\ldots,{\psi}_N(x,y)\}, we write .. _Eq:fem:sys:wT:ex:sum: .. math:: \tag{7} w = \sum_{j=0}^N c^{(w)}_j {\psi}_j,\quad T = T_0 + \sum_{j=0}^N c^{(T)}_j {\psi}_j{\thinspace .} Note that :math:`w` and :math:`T` in :ref:`(5) `-:ref:`(6) ` denote the exact solution of the PDEs, while :math:`w` and :math:`T` :ref:`(7) ` are the discrete functions that approximate the exact solution. It should be clear from the context whether a symbol means the exact or approximate solution, but when we need both at the same time, we use a subscript e to denote the exact solution. Variational form of each individual PDE --------------------------------------- Inserting the expansions :ref:`(7) ` in the governing PDEs, results in a residual in each equation, .. _Eq:fem:sys:wT:ex:weq:R: .. math:: \tag{8} R_w = \mu \nabla^2 w + \beta, .. _Eq:fem:sys:wT:ex:Teq:R: .. math:: \tag{9} R_T = \kappa\nabla^2 T + \mu ||\nabla w||^2 {\thinspace .} A Galerkin method demands :math:`R_w` and :math:`R_T` do be orthogonal to :math:`V`: .. math:: \int_\Omega R_w v {\, \mathrm{d}x} &=0\quad\forall v\in V,\\ \int_\Omega R_T v {\, \mathrm{d}x} &=0\quad\forall v\in V {\thinspace .} Because of the Dirichlet conditions, :math:`v=0` on :math:`\partial\Omega`. We integrate the Laplace terms by parts and note that the boundary terms vanish since :math:`v=0` on :math:`\partial\Omega`: .. _Eq:fem:sys:wT:ex:w:vf1: .. math:: \tag{10} \int_\Omega \mu \nabla w\cdot\nabla v {\, \mathrm{d}x} = \int_\Omega \beta v{\, \mathrm{d}x} \quad\forall v\in V, .. _Eq:fem:sys:wT:ex:T:vf1: .. math:: \tag{11} \int_\Omega \kappa \nabla T\cdot\nabla v {\, \mathrm{d}x} = \int_\Omega \mu \nabla w\cdot\nabla w\, v{\, \mathrm{d}x} \quad\forall v\in V {\thinspace .} Compound scalar variational form -------------------------------- The alternative way of deriving the variational from is to introduce a test vector function :math:`\boldsymbol{v}\in\boldsymbol{V} = V\times V` and take the inner product of :math:`\boldsymbol{v}` and the residuals, integrated over the domain: .. math:: \int_{\Omega} (R_w, R_T)\cdot\boldsymbol{v} {\, \mathrm{d}x} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V} {\thinspace .} With :math:`\boldsymbol{v} = (v_0,v_1)` we get .. math:: \int_{\Omega} (R_w v_0 + R_T v_1) {\, \mathrm{d}x} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V} {\thinspace .} Integrating the Laplace terms by parts results in .. _Eq:fem:sys:wT:ex:wT:vf2: .. math:: \tag{12} \int_\Omega (\mu\nabla w\cdot\nabla v_0 + \kappa\nabla T\cdot\nabla v_1){\, \mathrm{d}x} = \int_\Omega (\beta v_0 + \mu\nabla w\cdot\nabla w\, v_1){\, \mathrm{d}x}, \quad\forall \boldsymbol{v}\in\boldsymbol{V} {\thinspace .} Choosing :math:`v_0=v` and :math:`v_1=0` gives the variational form :ref:`(10) `, while :math:`v_0=0` and :math:`v_1=v` gives :ref:`(11) `. With the inner product notation, :math:`(p,q) = \int_\Omega pq{\, \mathrm{d}x}`, we can alternatively write :ref:`(10) ` and :ref:`(11) ` as .. math:: (\mu\nabla w,\nabla v) &= (\beta, v) \quad\forall v\in V,\\ (\kappa \nabla T,\nabla v) &= (\mu\nabla w\cdot\nabla w, v)\quad\forall v\in V, or since :math:`\mu` and :math:`\kappa` are considered constant, .. _Eq:fem:sys:wT:ex:w:vf1i: .. math:: \tag{13} \mu (\nabla w,\nabla v) = (\beta, v) \quad\forall v\in V, .. _Eq:fem:sys:wT:ex:T:vf1i: .. math:: \tag{14} \kappa(\nabla T,\nabla v) = \mu(\nabla w\cdot\nabla w, v)\quad\forall v\in V {\thinspace .} Decoupled linear systems ------------------------ The linear systems governing the coefficients :math:`c_j^{(w)}` and :math:`c_j^{(T)}`, :math:`j=0,\ldots,N`, are derived by inserting the expansions :ref:`(7) ` in :ref:`(10) ` and :ref:`(11) `, and choosing :math:`v={\psi}_i` for :math:`i=0,\ldots,N`. The result becomes .. _Eq:fem:sys:wT:ex:linsys:w1: .. math:: \tag{15} \sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N, .. _Eq:fem:sys:wT:ex:linsys:T1: .. math:: \tag{16} \sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N, .. _Eq:_auto2: .. math:: \tag{17} A^{(w)}_{i,j} = \mu(\nabla {\psi}_j,\nabla {\psi}_i), .. _Eq:_auto3: .. math:: \tag{18} b_i^{(w)} = (\beta, {\psi}_i), .. _Eq:_auto4: .. math:: \tag{19} A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j,\nabla {\psi}_i), .. _Eq:_auto5: .. math:: \tag{20} b_i^{(T)} = \mu((\sum_j c^{(w)}_j\nabla{\psi}_j)\cdot (\sum_k c^{(w)}_k\nabla{\psi}_k), {\psi}_i) {\thinspace .} It can also be instructive to write the linear systems using matrices and vectors. Define :math:`K` as the matrix corresponding to the Laplace operator :math:`\nabla^2`. That is, :math:`K_{i,j} = (\nabla {\psi}_j,\nabla {\psi}_i)`. Let us introduce the vectors .. math:: b^{(w)} &= (b_0^{(w)},\ldots,b_{N}^{(w)}),\\ b^{(T)} &= (b_0^{(T)},\ldots,b_{N}^{(T)}),\\ c^{(w)} &= (c_0^{(w)},\ldots,c_{N}^{(w)}),\\ c^{(T)} &= (c_0^{(T)},\ldots,c_{N}^{(T)}){\thinspace .} The system :ref:`(15) `-:ref:`(16) ` can now be expressed in matrix-vector form as .. _Eq:_auto6: .. math:: \tag{21} \mu K c^{(w)} = b^{(w)}, .. _Eq:_auto7: .. math:: \tag{22} \kappa K c^{(T)} = b^{(T)}{\thinspace .} We can solve the first system for :math:`c^{(w)}`, and then the right-hand side :math:`b^{(T)}` is known such that we can solve the second system for :math:`c^{(T)}`. Coupled linear systems ---------------------- Despite the fact that :math:`w` can be computed first, without knowing :math:`T`, we shall now pretend that :math:`w` and :math:`T` enter a two-way coupling such that we need to derive the algebraic equations as *one system* for all the unknowns :math:`c_j^{(w)}` and :math:`c_j^{(T)}`, :math:`j=0,\ldots,N`. This system is nonlinear in :math:`c_j^{(w)}` because of the :math:`\nabla w\cdot\nabla w` product. To remove this nonlinearity, imagine that we introduce an iteration method where we replace :math:`\nabla w\cdot\nabla w` by :math:`\nabla w_{-}\cdot\nabla w`, :math:`w_{-}` being the :math:`w` computed in the previous iteration. Then the term :math:`\nabla w_{-}\cdot\nabla w` is linear in :math:`w` since :math:`w_{-}` is known. The total linear system becomes .. _Eq:fem:sys:wT:ex:linsys:w2: .. math:: \tag{23} \sum_{j=0}^N A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(w,T)}_{i,j} c^{(T)}_j = b_i^{(w)},\quad i=0,\ldots,N, .. _Eq:fem:sys:wT:ex:linsys:T2: .. math:: \tag{24} \sum_{j=0}^N A^{(T,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(T,T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N, .. _Eq:_auto8: .. math:: \tag{25} A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j,{\psi}_i), .. _Eq:_auto9: .. math:: \tag{26} A^{(w,T)}_{i,j} = 0, .. _Eq:_auto10: .. math:: \tag{27} b_i^{(w)} = (\beta, {\psi}_i), .. _Eq:_auto11: .. math:: \tag{28} A^{(w,T)}_{i,j} = \mu((\nabla{\psi} w_{-})\cdot\nabla{\psi}_j), {\psi}_i), .. _Eq:_auto12: .. math:: \tag{29} A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j,{\psi}_i), .. _Eq:_auto13: .. math:: \tag{30} b_i^{(T)} = 0 {\thinspace .} This system can alternatively be written in matrix-vector form as .. _Eq:_auto14: .. math:: \tag{31} \mu K c^{(w)} = b^{(w)}, .. _Eq:_auto15: .. math:: \tag{32} L c^{(w)} + \kappa K c^{(T)} =0, with :math:`L` as the matrix from the :math:`\nabla w_{-}\cdot\nabla` operator: :math:`L_{i,j} = A^{(w,T)}_{i,j}`. The matrix-vector equations are often conveniently written in block form: .. math:: \left(\begin{array}{cc} \mu K & 0\\ L & \kappa K \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right), Note that in the general case where all unknowns enter all equations, we have to solve the compound system :ref:`(23) `-:ref:`(24) ` since then we cannot utilize the special property that :ref:`(15) ` does not involve :math:`T` and can be solved first. When the viscosity depends on the temperature, the :math:`\mu\nabla^2w` term must be replaced by :math:`\nabla\cdot (\mu(T)\nabla w)`, and then :math:`T` enters the equation for :math:`w`. Now we have a two-way coupling since both equations contain :math:`w` and :math:`T` and therefore must be solved simultaneously Th equation :math:`\nabla\cdot (\mu(T)\nabla w)=-\beta` is nonlinear, and if some iteration procedure is invoked, where we use a previously computed :math:`T_{-}` in the viscosity (:math:`\mu(T_{-})`), the coefficient is known, and the equation involves only one unknown, :math:`w`. In that case we are back to the one-way coupled set of PDEs. We may also formulate our PDE system as a vector equation. To this end, we introduce the vector of unknowns :math:`\boldsymbol{u} = (u^{(0)},u^{(1)})`, where :math:`u^{(0)}=w` and :math:`u^{(1)}=T`. We then have .. math:: \nabla^2 \boldsymbol{u} = \left(\begin{array}{cc} -{\mu}^{-1}{\beta}\\ -{\kappa}^{-1}\mu \nabla u^{(0)}\cdot\nabla u^{(0)} \end{array}\right){\thinspace .} Different function spaces for the unknowns ========================================== .. index:: mixed finite elements It is easy to generalize the previous formulation to the case where :math:`w\in V^{(w)}` and :math:`T\in V^{(T)}`, where :math:`V^{(w)}` and :math:`V^{(T)}` can be different spaces with different numbers of degrees of freedom. For example, we may use quadratic basis functions for :math:`w` and linear for :math:`T`. Approximation of the unknowns by different finite element spaces is known as *mixed finite element methods*. We write .. math:: V^{(w)} &= \hbox{span}\{{\psi}_0^{(w)},\ldots,{\psi}_{N_w}^{(w)}\},\\ V^{(T)} &= \hbox{span}\{{\psi}_0^{(T)},\ldots,{\psi}_{N_T}^{(T)}\} {\thinspace .} The next step is to multiply :ref:`(5) ` by a test function :math:`v^{(w)}\in V^{(w)}` and :ref:`(6) ` by a :math:`v^{(T)}\in V^{(T)}`, integrate by parts and arrive at .. _Eq:fem:sys:wT:ex:w:vf3: .. math:: \tag{33} \int_\Omega \mu \nabla w\cdot\nabla v^{(w)} {\, \mathrm{d}x} = \int_\Omega \beta v^{(w)}{\, \mathrm{d}x} \quad\forall v^{(w)}\in V^{(w)}, .. _Eq:fem:sys:wT:ex:T:vf3: .. math:: \tag{34} \int_\Omega \kappa \nabla T\cdot\nabla v^{(T)} {\, \mathrm{d}x} = \int_\Omega \mu \nabla w\cdot\nabla w\, v^{(T)}{\, \mathrm{d}x} \quad\forall v^{(T)}\in V^{(T)} {\thinspace .} The compound scalar variational formulation applies a test vector function :math:`\boldsymbol{v} = (v^{(w)}, v^{(T)})` and reads .. _Eq:fem:sys:wT:ex:wT:vf3: .. math:: \tag{35} \int_\Omega (\mu\nabla w\cdot\nabla v^{(w)} + \kappa\nabla T\cdot\nabla v^{(T)}){\, \mathrm{d}x} = \int_\Omega (\beta v^{(w)} + \mu\nabla w\cdot\nabla w\, v^{(T)}){\, \mathrm{d}x}, valid :math:`\forall \boldsymbol{v}\in\boldsymbol{V} = V^{(w)}\times V^{(T)}`. The associated linear system is similar to :ref:`(15) `-:ref:`(16) ` or :ref:`(23) `-:ref:`(24) `, except that we need to distinguish between :math:`{\psi}_i^{(w)}` and :math:`{\psi}_i^{(T)}`, and the range in the sums over :math:`j` must match the number of degrees of freedom in the spaces :math:`V^{(w)}` and :math:`V^{(T)}`. The formulas become .. _Eq:fem:sys:wT:ex:linsys:w1:mixed: .. math:: \tag{36} \sum_{j=0}^{N_w} A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N_w, .. _Eq:fem:sys:wT:ex:linsys:T1:mixed: .. math:: \tag{37} \sum_{j=0}^{N_T} A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N_T, .. _Eq:_auto16: .. math:: \tag{38} A^{(w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},{\psi}_i^{(w)}), .. _Eq:_auto17: .. math:: \tag{39} b_i^{(w)} = (\beta, {\psi}_i^{(w)}), .. _Eq:_auto18: .. math:: \tag{40} A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},{\psi}_i^{(T)}), .. _Eq:_auto19: .. math:: \tag{41} b_i^{(T)} = \mu(\nabla w_{-}, {\psi}_i^{(T)}) {\thinspace .} In the case we formulate one compound linear system involving both :math:`c^{(w)}_j`, :math:`j=0,\ldots,N_w`, and :math:`c^{(T)}_j`, :math:`j=0,\ldots,N_T`, :ref:`(23) `-:ref:`(24) ` becomes .. _Eq:fem:sys:wT:ex:linsys:w2b: .. math:: \tag{42} \sum_{j=0}^{N_w} A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^{N_T} A^{(w,T)}_{i,j} c^{(T)}_j = b_i^{(w)},\quad i=0,\ldots,N_w, .. _Eq:fem:sys:wT:ex:linsys:T2b: .. math:: \tag{43} \sum_{j=0}^{N_w} A^{(T,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^{N_T} A^{(T,T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N_T, .. _Eq:_auto20: .. math:: \tag{44} A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},{\psi}_i^{(w)}), .. _Eq:_auto21: .. math:: \tag{45} A^{(w,T)}_{i,j} = 0, .. _Eq:_auto22: .. math:: \tag{46} b_i^{(w)} = (\beta, {\psi}_i^{(w)}), .. _Eq:_auto23: .. math:: \tag{47} A^{(w,T)}_{i,j} = \mu (\nabla w_{-}\cdot\nabla{\psi}_j^{(w)}), {\psi}_i^{(T)}), .. _Eq:_auto24: .. math:: \tag{48} A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},{\psi}_i^{(T)}), .. _Eq:_auto25: .. math:: \tag{49} b_i^{(T)} = 0 {\thinspace .} The corresponding block form .. math:: \left(\begin{array}{cc} \mu K^{(w)} & 0\\ L & \kappa K^{(T)} \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right), has square and rectangular block matrices: :math:`K^{(w)}` is :math:`N_w\times N_w`, :math:`K^{(T)}` is :math:`N_T\times N_T`, while :math:`L` is :math:`N_T\times N_w`, Computations in 1D ================== .. 2DO .. show analytical solution in [0,H] .. use global polynomials (x^i(H-x)), exact sol .. compute uncoupled and coupled discrete systems, N=4 .. note: coupled can use exact w_{-} .. P1-P1, n elements, 2 elements as special case .. uncoupled and coupled (can use exercises for variants) .. P1 for w, P2 for T or P4 for T .. similar computations for circle can be done as project .. extensions to time-dep versions in projects .. any geophysical applications? flowing ice sheet ("half channel") and .. temp gradient through, check with Jed .. the time-dependent system can be introduced in diffusion and .. a finite difference scheme can be devised We can reduce the system :ref:`(5) `-:ref:`(6) ` to one space dimension, which corresponds to flow in a channel between two flat plates. Alternatively, one may consider flow in a circular pipe, introduce cylindrical coordinates, and utilize the radial symmetry to reduce the equations to a one-dimensional problem in the radial coordinate. The former model becomes .. _Eq:fem:sys:wT:ex1D:weq: .. math:: \tag{50} \mu w_{xx} = -\beta, .. _Eq:fem:sys:wT:ex1D:Teq: .. math:: \tag{51} \kappa T_{xx} = - \mu w_x^2, while the model in the radial coordinate :math:`r` reads .. _Eq:fem:sys:wT:ex1Dr:weq: .. math:: \tag{52} \mu\frac{1}{r}\frac{d}{dr}\left( r\frac{dw}{dr}\right) = -\beta, .. _Eq:fem:sys:wT:ex1Dr:Teq: .. math:: \tag{53} \kappa \frac{1}{r}\frac{d}{dr}\left( r\frac{dT}{dr}\right) = - \mu \left( \frac{dw}{dr}\right)^2 {\thinspace .} The domain for :ref:`(50) `-:ref:`(51) ` is :math:`\Omega = [0,H]`, with boundary conditions :math:`w(0)=w(H)=0` and :math:`T(0)=T(H)=T_0`. For :ref:`(52) `-:ref:`(53) ` the domain is :math:`[0,R]` (:math:`R` being the radius of the pipe) and the boundary conditions are :math:`du/dr = dT/dr =0` for :math:`r=0`, :math:`u(R)=0`, and :math:`T(R)=T_0`. .. 2DO **Calculations to be continued...**