I am trying to solve a 1D initial boundary value problem in MATLAB using Finite Elements with time stepping, for the purpose of learning scientific computing and to build up to more difficult problems. The initial boundary value problem is:
Find $u(x,t)$ such that: \begin{align*} u_{xxxx} + u_{tt} = f(x,t), \\ u(0,t) = 0, \quad &u(1,t) = 0, \\ u_x(0,t) = \frac{\pi t^2}{100}, \quad &u_x(1,t) = -\frac{\pi t^2}{100}, \\ u(x,0) = 0, \quad & u(1,t) = 0, \\ f(x,t) = 0.01\sin(\pi x)(\pi^4 t^2 + 2), \end{align*} for $x \in \Omega = [0,1]$ and $t \in \Omega_t = [0,t_{end}]$.
The solution to this problem is $u(x,t) = 0.01 t^2 \sin(\pi x)$. The weak/variational form is:
Find $u \in S$ such that for all $v \in T$: \begin{align*} \int_0^1 u_{xx}v_{xx} dx + \int_0^1 u_{tt} v dx= \int_0^1 v f dx, \end{align*} where $$T = \{ v: \int |v|^2 < \infty, \int|v_x|^2 < \infty, \int |v_{xx}|^2 < \infty , v(0) =v(1)=v_x(0)=v_x(1) = 0 \},$$ and $$ S = \{u:\int |u|^2 < \infty, \int|u_x|^2 < \infty, \int |u_{xx}|^2 < \infty , u(0,t) = u(1,t) = 0, u_x(0,t) = -u_x(1,t) = \frac{\pi t^2}{100} \}.$$
This weak form will lead to the final finite-element system: $$ A_{22} u(t) + A_{00} u_{tt}(t) = f_0(t) + b_{A_{22}}(t) + b_{A_{00}(t)}. $$ Here, the entries of the matrix $A_{22}$ can be computed "locally" from $\int_0^1 u_{xx} v_{xx} dx$ and the entries of the matrix $A_{00}$ can be computed locally from $\int_0^1 u_{tt} v dx$, then the locally computed entries can be assembled into the final "global" matrices. The entries of the time-dependent vector $f_0 (t)$ can be computed locally from $\int_0^1 vf dx$, and can once again be assembled into a final global "load vector". Also the time-dependent vectors $b_{A_{22}}(t)$ and $b_{A_{00}}(t)$ come from enforcing the essential boundary conditions.
I have already solved a very similar steady-state problem by introducing a Mesh for the domain, making relevant data structures and functions to access $C^1$ finite element basis functions, integrating the products of these basis functions to compute the relevant matrices/vectors, and solving the system of equations. Doing this resulted in an array of coefficients of $u(x)$ and $u_x(x)$ at the Mesh points, so that the solution to the variational problem could be approximated using the $C^1$ finite element basis.
Problem: Solve the time-dependent problem presented here using a fully discrete Finite Element scheme.
My attempt:
For the fully discrete problem, we introduce a partition for the time interval $[0, t_{end}]$: $$ 0 = t_0 < t_1 < ... < t_m = t_{end}, $$ with time step sizes $\tau_{k} = (t_k - t_{k-1})$. Then our goal is to find an $m$ number of vectors $u^k \approx u(t_k)$, i.e., each vector contains the coefficients of $u(x,t_k)$ and $u_x(x,t_k)$ at the $k$-th time step (so that the solution at each time step can be interpolated using the basis functions). I specifically wanted to use the Crank-Nicolson method. First I rewrite my second-order system as a first-order system by introducing variables $y_1(t) = u(t)$ and $y_2(t) = u_t (t)$, then differentiating these variables gives us: \begin{align*} \dot{y}(t) = \begin{bmatrix} \dot{y}_1(t) \\ \dot{y}_2(t) \end{bmatrix} = \begin{bmatrix} y_2(t) \\ A_{00}^{-1} (-A_{22} y_1(t) + f_0(t) + b_{A_{22}}(t) + b_{A_{00}}(t)) \end{bmatrix}, \end{align*} which is a first order system. Letting $t_{k - 1/2} = (t_{k - 1} + t_k)/2$, we use the following approximations for $y$ and $y'$ at the $t_{k - 1/2}$-th time step: \begin{align*} y(t_{k - 1/2}) = \begin{bmatrix} y_1 (t_{k - 1/2}) \\ y_2 (t_{k - 1/2}) \end{bmatrix} \approx \frac{y(t_{k-1}) + y(t_k)}{2}, \\ \dot{y}(t_{k - 1/2}) \approx \frac{y(t_k) - y(t_{k-1})}{\tau_k}. \end{align*} Substituting these into the first-order system and rearranging to obtain the $y_1 (t_k)$ and $y_2 (t_k)$ on the right hand side gives \begin{align*} \begin{bmatrix} y_1 (t_k)/\tau_k - y_2 (t_k)/2 \\ A_{00} y_2(t_k)/\tau_k + A_{22} y_1 (t_k)/2 \end{bmatrix} = \begin{bmatrix} y_1(t_{k-1})/\tau_k + y_2(t_{k - 1})/2 \\ A_{00} y_2 (t_{k-1})/\tau_k - A_{22} y_1 (t_{k-1})/2 + f_0(t_{k - 1/2}) + b_{A_{22}}(t_{k - 1/2}) + b_{A_{00}}(t_{k - 1/2}) \end{bmatrix}. \end{align*} Is this the correct formula I should use, or is there anything I am missing? The initial conditions here are just $y_1(0) = 0$ and $y_2(0) = 0$, and I could just form $f_0 (t_k)$ at each time step, so it seems like things should run smoothly. For example, for the first time step $k = 1$ I would have \begin{align*} \begin{bmatrix} y_1 (t_1)/\tau_1 - y_2 (t_1)/2 \\ A_{00} y_2(t_1)/\tau_1 + A_{22} y_1 (t_1)/2 \end{bmatrix} = \begin{bmatrix} 0 \\ f_0(t_{1/2}) + b_{A_{22}}(t_{1/2}) + b_{A_{00}}(t_{1/2}) \end{bmatrix}, \end{align*} and here the only unknowns are $y_1(t_1)$ and $y_2(t_2)$.
EDIT: It turns out that there was a small mistake in my original MATLAB code and everything seems to be working fine. Should I delete this question, or what?
