# Method of lines for inhomogeneous Dirichlet conditions

I understand how to set up the boundary conditions for a steady state problem discretized by Galerkin method, for a time dependent PDE below, $$\frac{\partial}{\partial t} u = c\nabla^2 u + a\nabla u + b,~~~ u_0= 1, u'_1=0$$

Supposing I have correctly set up the Dirichlet and Neumann conditions for the steady state part $c\nabla^2 u + a\nabla u + b=0$, which leads to an algebraic system $C u + A u + B=0$.

Then turning to the time dependent problem, I have an additional "mass" matrix at left hand side, $$Mu=\frac{\partial}{\partial t} u, ~~~~~~~~ M=M(v,u)$$

The problem I encountered is: If I don't fix the Dirichlet end at each time iteration(force it to be the constant, seems a little "dirty"), the solution would blow up at the Dirichlet boundary in the end. also some oscillation can be observed(depends on the specific $\theta$ chosen) near the Dirichlet boundary during the integration.

My questions are:

1) I don't have to modify the matrix $C, A, B$ (boundary part), right?

2) what is the Boundary condition for $M$ (if the $\theta$ method is to be implemented), also Dirichlet and Neumann respectively? I.e. use the same basis function as the time-indepedent part.

• @ChristianClason Thanks! I want to implement the $\theta$ method, in particular, does the specific value of $\theta$ affect the boundary condition? – lorniper Jun 28 '14 at 20:27
• From your comments below, it seems that your problem is specifically how to treat the inhomogeneous Dirichlet conditions. I edited the title of your question to make this clear; I hope you don't mind. (Otherwise feel free to roll back my edit.) – Christian Clason Jun 29 '14 at 10:17

In a Galerkin approach for stationary PDEs, you write your unknown solution as a linear combination of "nice" basis functions $\{\phi_j\}_{j=1,\dots,n}$, $$u(x) \approx \sum_{j=1}^n u_j \phi_j(x),$$ and project the equation onto the space spanned by (these or other) functions, exploiting the linearity of the differential operators and the inner product. In your case, you get $$\sum_{j=1}^n u_j (c\nabla^2 \phi_j(x)+a\nabla \phi_j(x) + b,\psi_i(x))=0\qquad\text{for all }1\leq i\leq n.$$ This yields a set of linear equations for the unknown coefficients $u_j$.

In the finite element method, boundary conditions are implemented differently for Dirichlet and for Neumann conditions. Homogeneous Dirichlet conditions are built into the choice of the basis: $\phi_j(0)=\psi_i(0)=0$ for all $1\leq i,j\leq n$. The Neumann conditions, on the other hand, become part of the operator: You integrate by parts in the first term of the sum to obtain (I am leaving out the lower order terms since they do not play such a role here) $$\sum_{j=1}^n u_j \left(\phi_j'(1)\psi_i(1)-(c\nabla \phi_j(x),\nabla\psi_i(x)) \right)\qquad\text{for all }1\leq i\leq n$$ (using that all boundary terms for $x=0$ vanish due to the choice of $\psi$). Of course, then you use that $$\sum_{j=1}^n u_j \phi_j'(1)\psi_i(1) = u'(1)\psi_i(1)$$ so that becomes part of your right-hand side. In your case, $u'(1)=0$ so the whole term vanishes.

In the method of lines, you approximate your unknown solution $u(t,x)$ as $$u(t,x)\approx \sum_{j=1}^n u_j(t) \phi_j(x),$$ i.e., the whole space dependence is in the basis functions, while the whole time dependence is in the scalar coefficients. Your projected equation is now $$\sum_{j=1}^n u_j(t) (c\nabla^2 \phi_j(x)+a\nabla \phi_j(x) + b,\psi_i(x)) = \sum_{j=1}^n \frac{d}{dt}u_j(t)(\phi_j(x),\psi_i(x))\qquad\text{for all }1\leq i\leq n.$$ Now you do the same steps; since the inner products on the left-hand side are exactly the same as before, you get exactly the same matrices. On the right-hand side, you have a term just like a standard mass matrix $M$ with entries $M_{ij}=(\phi_j(x),\psi_i(x))$, where you respect the Dirichlet conditions by the choice of basis functions, but you don't have to worry about Neumann conditions. Note that the time dependence is purely in the coefficients, so for the (spatial) mass matrix it makes no difference whether this comes from approximating $u$ or $\partial_t u$.

The way inhomogeneous Dirichlet conditions are treated mathematically is to reduce the problem to the homogeneous case: You write $$u(t,x) = g(t,x) + \tilde u(t,x)$$ for some fixed $g$ with $g(t,0)=1$ and $g'(t,1)=0$. Then, $\tilde u$ satisfies $$c\nabla^2 \tilde u + a\nabla \tilde u = -b -c\nabla^2 g -a\nabla g,\qquad \tilde u(0) = \tilde u'(1) = 0,$$ and you can proceed as above to compute $\tilde u(t,x)$ and then $u(t,x) = \tilde u(t,x)+g(t,x)$. If you work with the standard hat functions, you could simply take $g(t,x) = \phi_1(x)$.

If the initial condition (which you have not stated) satisfies $u_0(0)=1$, you could also replace the first equation of the system of ODEs by $\partial_t u_1 = 0$, i.e., replace the first row in the mass matrix by $(1,0,\dots,0)$, the first row of the remaining matrices by $(0,\dots,0)$, and the first component of the right-hand side by $0$.

If memory serves correctly, the method of lines is discussed in Ern, Guermond: Theory and Practice of Finite Elements, Springer, 2004, Chapter 6.1.

• in the third equation, after you do integration by part, its $\psi(1)\nabla\phi(1)$ right? and why did you neglect the Neumann part for the convection term, doesn't something come out of the integration by part as well? – lorniper Jun 28 '14 at 21:42
• You usually only integrate by parts in the second-order term because this way you only need to have first derivatives of $u$ in the equation. Since you can't get rid of all derivatives easily, you might as well keep the convection term as is. – Christian Clason Jun 28 '14 at 23:26
• thanks! but still one problem, since the mass matrix is assembled by the chosen basis functions, on the Dirichlet end of "M", practically what do we do? still eliminate first line, and put "1" in the first diagonal entry? or some other tricks? – lorniper Jun 29 '14 at 8:04
• being back after a long time,sigh... now I wonder why I "don't have to worry about Neumann conditions.", is it because we have an additional degree of freedom for time, so the Neumann condition is already relaxed? – lorniper Dec 1 '14 at 20:39