I have a time-dependent two dimensional pde that I am solving by FEM. I have been successful solving it with forward Euler. Now, I am trying to solve with higher order methods such as RK4. It is fairly straight forward but the hard part is implementing the Dirichlet boundary condition. For example, with forward Euler $A \frac{dx}{dt}=b$ then with forward Euler $A \frac{x_{t+1}-x_t}{\Delta t } = b $. Thus, $A x_{t+1} = A x_t +\Delta t b$ which then the Dirichlet boundary condition can be implemented. But I am not sure how this works with RK4.
Could anyone give me an explanation how this is done or recommend me a good reading for this?
You might already have the answer but let me comment nevertheless. You can use the same approach as when deriving a formula for non-homogeneous Dirichlet boundary conditions. Let me give you an example.
Consider the heat equation $$ u_t=\Delta u. $$ Let $\phi_i,~i\in\{1,\cdots,N\}$ be your basis. Discretizing the weak formulation gives $$ M \dot{x} + A x = 0, $$ where $M_{ij} = (\phi_i,\phi_j)$, $A_{ij}=(\nabla \phi_i, \nabla \phi_j)$ and $x \in \mathbb{R}^N$ is the vector of unknown degrees-of-freedom. Next split your DOF vector in such a way that the DOF's corresponding to the Dirichlet boundary come last, i.e. $$ x=\begin{bmatrix} x_I \\ x_B \end{bmatrix}, $$ where $x_I$ consists of the interior (or free) DOFs and $x_B$ consists of the boundary DOFs (for which you should know the value already). This kind of splitting induces a similar splitting for the discrete weak formulation: $$ \begin{bmatrix} M_{II} & M_{IB} \\ M_{BI} & M_{BB} \end{bmatrix} \begin{bmatrix} \dot{x_I} \\ \dot{x_B} \end{bmatrix}+\begin{bmatrix} A_{II} & A_{IB} \\ A_{BI} & A_{BB} \end{bmatrix} \begin{bmatrix} x_I \\ x_B \end{bmatrix} = 0. $$ The first set of equations reads now $$ M_{II} \dot{x_I} + A_{II} x_I = -A_{IB} x_B - M_{IB} \dot{x_B}, $$ which you can solve for $x_I$ using any difference method. Note that the right hand side consists only of terms that you know a priori. The second term on the right hand side might be zero if your Dirichlet boundary condition does not depend on time.