I am reading a textbook on PDE (Introduction to Elliptic and Parabolic Partial Differential Equations by Z. Wu, J. Yin and C. Wang, written in Chinese) where a proof of the existence of weak solutions of heat equation by Rothe's method is given. Now I am trying to apply this method to the more general case: \begin{align} &\frac{\partial u}{\partial t}-\sum_{i,j=1}^n \frac{\partial}{\partial x_j} \left(a_{ij} \frac{\partial u}{\partial x_i}\right) + cu=f, \quad (x,t) \in Q_T, \quad (1) \\ &u(x,t)=0, \quad (x,t) \in \partial \Omega \times (0,T), \qquad (2) \\ &u(x,0)=u_0(x), \quad x \in \Omega, \qquad (3) \end{align} where $\Omega$ is a bounded open subset of $\mathbb{R}^n$ with a Lipschitz boundary, $T$ is a positive constant, $Q_T:=\Omega \times (0,T)$, $f \in L^2(Q_T)$, $u_0 \in H_0^1(\Omega)$, and there exist two positive constants $\lambda$ and $M$, such that \begin{gather} a_{ij}(x,t) \xi_i \xi_j \geq \lambda |\xi|^2, \quad \forall (x,t) \in Q_T, \ \forall \xi=(\xi_1,\dots,\xi_n) \in \mathbb{R}^n, \\ \sum_{i,j=1}^n \|a_{ij}\|_{L^\infty(Q_T)}+\|c\|_{L^\infty(Q_T)} \leq M. \end{gather} Denote $\overset{\circ}{C^\infty}(\overline{Q_T})$ as the set of smooth funtions which vanish near the lateral boundary $\partial \Omega \times (0,T)$ and $\ \ \overset{\circ}{\vphantom{W}}\mspace{-12mu}W^{k_1,k_2}_p(Q_T)$ the closure of $\overset{\circ}{C^\infty}(\overline{Q_T})$ in the Sobolev space $W^{k_1,k_2}_p(Q_T)$. We say $u \in \ \ \overset{\circ}{\vphantom{W}}\mspace{-12mu}W^{1,1}_2(Q_T)$ is a weak solution of (1)-(3), if for any $\varphi \in \ \ \overset{\circ}{\vphantom{W}}\mspace{-12mu}W^{1,0}_2(Q_T)$ there hold \begin{equation} \iint_{Q_T} \left(\frac{\partial u}{\partial t} \varphi + \sum_{i,j=1}^n a_{ij} \frac{\partial u}{\partial x_i} \frac{\partial \varphi}{\partial x_j} + cu \varphi \right) \, \mathrm{d} x \mathrm{d} t = \iint_{Q_T} f \varphi \, \mathrm{d} x \mathrm{d} t \end{equation} and \begin{equation} \gamma u(x,0)=u_0(x), \quad a.e. x \in \Omega, \end{equation} where $\gamma u(x,0)$ means the trace of $u$ on the bottom boundary $\Omega \times \{t=0\}$.
My proof is as follows, which extends the original proof step by step:
Step 1: Discretize the time variable $t$ and construct the approximate solutions $\{u^m\}$.
For each positive integer $m$ and function $w(x,t)$, denote $w^{m,l}(x):=w(x,lh)$ $(l=0,1,\dots,m)$, where $h=T/m$. For each $l$, consider the approximate equation:
\begin{equation}
\frac{u^{m,l}-u^{m,l-1}}{h}-\sum_{i,j=1}^n \frac{\partial}{\partial x_j} \left(a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i}\right) + c^{m,l} u^{m,l} = f^{m,l}, \quad x \in \Omega,
\end{equation}
or equivalently,
\begin{equation}
-\sum_{i,j=1}^n \frac{\partial}{\partial x_j} \left(a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i}\right) + \left(\frac{1}{h}+c^{m,l}\right) u^{m,l} = f^{m,l}+\frac{u^{m,l-1}}{h}, \quad x \in \Omega. \quad (4)
\end{equation}
Suppose $h$ is sufficiently small, so that $1/h+c^{m,l} \geq 0$ always holds, then being aware of the fact that $u^{m,0}=u_0 \in H_0^1(\Omega)$ and by using the standard theory for elliptic PDEs and mathematical induction, we deduce that there exist $u^{m,1},u^{m,2},\dots \in H_0^1(\Omega)$ which solve (4) for $l=1,2,\dots$ respectively, that is, for any $\varphi \in H_0^1(\Omega)$, we have
\begin{equation}
\int_\Omega \left(\frac{u^{m,l}-u^{m,l-1}}{h} \varphi + \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial \varphi}{\partial x_j} + c^{m,l} u^{m,l} \varphi \right) \, \mathrm{d} x = \int_\Omega f^{m,l} \varphi \, \mathrm{d} x, \quad l=1,\dots,m. \quad (5)
\end{equation}
Now let
\begin{gather}
u^m(x,t):=\sum_{l=1}^m \chi_{m,l}(t) u^{m,l-1}(x) + \sum_{l=1}^m \chi_{m,l}(t) \eta_{m,l}(t) (u^{m,l}(x)-u^{m,l-1}(x)), \\
w^m(x,t):=\sum_{l=1}^m \chi_{m,l}(t) u^{m,l}(x), \quad f^m(x,t):=\sum_{l=1}^m \chi_{m,l}(t) f^{m,l}(x),
\end{gather}
where $\chi_{m,l}(t)$ is the indicator function of the interval $[(l-1)h,lh)$, and
\begin{equation}
\eta_{m,l}(t):=
\begin{cases}
\frac{t}{h}-(l-1), \quad t \in [(l-1)h,lh), \\
0, \quad \text{elsewhere}.
\end{cases}
\end{equation}
The graphs of $u^m,w^m,f^m$ are basically as follow:
for each fixed $x \in \Omega$ (the position on the vertical axis doesn't mean which one is greater or smaller)
Then by (5), for each $\varphi \in H_0^1(\Omega)$, we have \begin{equation} \int_\Omega \left(\frac{\partial u^m}{\partial t} + \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial w^m}{\partial x_i} \frac{\partial \varphi}{\partial x_j} + c^{m,l} w^m \varphi \right) \, \mathrm{d} x = \int_\Omega f^m \varphi \, \mathrm{d} x, \quad t \in [0,T]. \end{equation}
Step 2: Estimate for $u^m$.
Set $\varphi=u^{m,l}-u^{m,l-1}$ in (5), then we have \begin{equation} \frac{1}{h} \left\|u^{m,l}-u^{m,l-1} \right\|^2_{L^2(\Omega)} = -\int_\Omega \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial}{\partial x_j} (u^{m,l}-u^{m,l-1}) \, \mathrm{d} x + \int_\Omega (f^{m,l}-c^{m,l} u^{m,l}) (u^{m,l}-u^{m,l-1}) \, \mathrm{d} x. \quad (6) \end{equation} We now scrutinize each term in (6). First we have \begin{equation} \begin{split} -\int_\Omega \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial}{\partial x_j} (u^{m,l}-u^{m,l-1}) \, \mathrm{d} x =& -\int_\Omega \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial u^{m,l}}{\partial x_j} \, \mathrm{d} x + \int_\Omega \sum_{i,j=1}^n a_{ij}^{m,l} \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial u^{m,l-1}}{\partial x_j} \, \mathrm{d} x \\ \leq& -\lambda \|\nabla u^{m,l}\|^2_{L^2(\Omega)} + \int_\Omega \sum_{i,j=1}^n M \frac{\partial u^{m,l}}{\partial x_i} \frac{\partial u^{m,l-1}}{\partial x_j} \, \mathrm{d} x \\ \leq& -\lambda \|\nabla u^{m,l}\|^2_{L^2(\Omega)} + \frac{n \varepsilon}{2} \|\nabla u^{m,l}\|^2_{L^2(\Omega)} + \frac{nM^2}{2 \varepsilon} \|\nabla u^{m,l-1}\|^2_{L^2(\Omega)}, \quad (7) \end{split} \end{equation} where $\varepsilon$ is a small positive number. Also we have \begin{equation} \begin{split} \int_\Omega (f^{m,l}-c^{m,l} u^{m,l}) (u^{m,l}-u^{m,l-1}) \, \mathrm{d} x \leq& \frac{h}{2} (\|f^{m,l}\|^2_{L^2(\Omega)}+M^2 \|u^{m,l}\|^2_{L^2(\Omega)}) + \frac{1}{2h} \|u^{m,l}-u^{m,l-1}\|^2_{L^2(\Omega)} \\ \leq& \frac{h}{2} \|f^{m,l}\|^2_{L^2(\Omega)} + \frac{\mu^2 M^2 h}{2} \|\nabla u^{m,l}\|^2_{L^2(\Omega)} + \frac{1}{2h} \|u^{m,l}-u^{m,l-1}\|^2_{L^2(\Omega)}, \quad (8) \end{split} \end{equation} where $\mu>0$ is the constant in the Poincare's inequality which depends only on $n$, that is, $\|u\|_{L^2(\Omega)} \leq \mu \|\nabla u\|_{L^2(\Omega)}$ $(\forall u \in H_0^1(\Omega))$.
Combining (6)-(8), we get \begin{equation} \frac{1}{h} \left\| u^{m,l}-u^{m,l-1} \right\|^2_{L^2(\Omega)} + (2 \lambda-\varepsilon n-\mu^2 M^2 h) \|\nabla u^{m,l}\|^2_{L^2(\Omega)} \leq \frac{nM^2}{\varepsilon} \|\nabla u^{m,l-1}\|^2_{L^2(\Omega)} + h \|f^{m,l}\|^2_{L^2(\Omega)}. \quad (9) \end{equation} Dropping the first term on the left-hand-side of (9), we get \begin{equation} (2 \lambda-\varepsilon n-\mu^2 M^2 h) \|\nabla u^{m,l}\|^2_{L^2(\Omega)} \leq \frac{nM^2}{\varepsilon} \|\nabla u^{m,l-1}\|^2_{L^2(\Omega)} + h \|f^{m,l}\|^2_{L^2(\Omega)}. \end{equation} Set $\varepsilon,h$ small enough, so that $2 \lambda-\varepsilon n-\mu^2 M^2 h>0$, and by dividing both sides by $2 \lambda-\varepsilon n-\mu^2 M^2 h$, we get \begin{equation} \|\nabla u^{m,l}\|^2_{L^2(\Omega)} \leq C (\|\nabla u^{m,l-1}\|^2_{L^2(\Omega)} + h \|f^{m,l}\|^2_{L^2(\Omega)}), \quad (10) \end{equation} where $C$ is a positive constant which depends on $n,\lambda,M,\varepsilon,h$.
Iterate (10) by $l$ times, and we get \begin{equation} \|\nabla u^{m,l}\|^2_{L^2(\Omega)} \leq C \left( \|\nabla u_0\|^2_{L^2(\Omega)} + h \sum_{k=1}^l \|f^{m,k}\|^2_{L^2(\Omega)} \right) \leq C \left( \|\nabla u_0\|^2_{L^2(\Omega)} + h \sum_{k=1}^m \|f^{m,k}\|^2_{L^2(\Omega)} \right) = C ( \|\nabla u_0\|^2_{L^2(\Omega)} + \|f^m\|^2_{L^2(Q_T)} ) = C M_m, \quad l=1,\dots,m, \quad (11) \end{equation} where $M_m:=\|\nabla u_0\|^2_{L^2(\Omega)} + \|f^m\|^2_{L^2(Q_T)}$, and $C$ is a positive constant which depends on $n,\lambda,M,\varepsilon,h,m$.
As a simple calculation shows that \begin{equation} \nabla u^m=\sum_{l=1}^m \chi_{m,l} (\nabla u^{m,l-1}+\eta_{m,l} (\nabla u^{m,l}-\nabla u^{m,l-1})), \end{equation} we derive from (11) that \begin{align} \left\|\nabla u^{m}\right\|_{L^{2}\left(Q_{T}\right)}^{2} &=\sum_{l=1}^{m} \int_{0}^{T} \chi_{m, l} \int_{\Omega}\left|\left(1-\eta_{m,l}\right) \nabla u^{m,l-1}+\eta_{m,l} \nabla u^{m, l}\right|^{2} \, \mathrm{d} x \mathrm{d} t \\ & \leq 2 \sum_{l=1}^{m} \int_{0}^{T} \chi_{m, l}\left(\left\|\nabla u^{m, l-1}\right\|_{L^{2}(\Omega)}^{2}+\left\|\nabla u^{m, l}\right\|_{L^{2}(\Omega)}^{2}\right) \, \mathrm{d} t \\ & \leq 2 \sum_{l=1}^{m} h\left(\left\|\nabla u^{m, l-1}\right\|_{L^{2}(\Omega)}^{2}+\left\|\nabla u^{m, l}\right\|_{L^{2}(\Omega)}^{2}\right) \leq 4 T C M_m, \quad (12) \end{align} where $C$ is a positive constant which depends on $n,\lambda,M,\varepsilon,h,m$.
The rest of the proof in the special case of heat equation (i.e. $a_{ij}=1$ for $i=j$, $a_{ij}=0$ for $i \neq j$, and $c=0$) is basically as follows: consider the special case where $f \in C(\overline{Q_T})$, so that $\|f^m-f\|_{L^2(Q_T)} \to 0$ as $m \to \infty$ (which can be proved by using the definition of Riemann integral), and therefore \begin{equation} \lim_{m \to \infty} M_m = \|\nabla u_0\|^2_{L^2(\Omega)}+\|f\|^2_{L^2(Q_T)}, \end{equation} which implies that $\{M_m\}_{m=1}^\infty$ is a bounded sequence, and so is $\{\|\nabla u^{m}\|^2_{L^2(Q_T)}\}_{m=1}^\infty$ according to (12). So here comes the first question, as in the special case of heat equation, the constant $C$ in (12) is independent of $m$ (actually $C=1$), this deduction is correct, however, in my proof, $C$ is dependent on $m$, therefore the boundedness of $\{M_m\}_{m=1}^\infty$ doesn't imply the boundedness of $\{\|\nabla u^{m}\|^2_{L^2(Q_T)}\}_{m=1}^\infty$. Next, in the original proof, (9), which reduces to \begin{equation} \frac{1}{h} \left\|u^{m,l}-u^{m,l-1}\right\|^2_{L^2(\Omega)} + \|\nabla u^{m,l}\|^2_{L^2(\Omega)} \leq \|\nabla u^{m,l-1}\|^2_{L^2(\Omega)} + h \|f^{m,l}\|^2_{L^2(\Omega)}, \quad (13) \end{equation} is summated to get \begin{equation} \frac{1}{h} \sum_{l=1}^m \left\|u^{m,l}-u^{m,l-1}\right\|^2_{L^2(\Omega)} \leq \|\nabla u_0\|^2_{L^2(\Omega)} + h \sum_{l=1}^m \|f^{m,l}\|^2_{L^2(\Omega)} = \|\nabla u_0\|^2_{L^2(\Omega)} + \|f^m\|^2_{L^2(Q_T)} = M_m, \quad (14) \end{equation} and as a simple calculation shows that \begin{equation} \frac{\partial u^m}{\partial t}=\frac{1}{h} \sum_{l=1}^m \chi_{m,l} (u^{m,l}-u^{m,l-1}), \end{equation} we get \begin{equation} \left\|\frac{\partial u^m}{\partial t}\right\|^2_{L^2(Q_T)} = \frac{1}{h^2} \sum_{l=1}^m h \left\|u^{m,l}-u^{m,l-1}\right\|^2_{L^2(\Omega)} \leq M_m, \end{equation} and therefore $\left\{\left\|\frac{\partial u^m}{\partial t}\right\|^2_{L^2(Q_T)}\right\}_{m=1}^\infty$ is bounded. In this way it is proved that $\{u^m\}$ is a bounded sequence in $W^{1,1}_2(Q_T)$, and so has a subsequence which converges weakly to a function, say $u$, in $W^{1,1}_2(Q_T)$. The rest is basically to prove that $u$ solves the problem (1)-(3) in weak sense and discuss the case where $f \in L^2(Q_T)$ rather than $f \in C(\overline{Q_T})$.
Then here comes the second question. The motivation of summating (13) to get (14) is to cancel out $\|\nabla u^{m,l}\|^2_{L^2(\Omega)}$ for $l=1,\dots,m$, but in my proof, the coefficient of $\|\nabla u^{m,l}\|^2_{L^2(\Omega)}$ and $\|\nabla u^{m,l-1}\|^2_{L^2(\Omega)}$ in (9) don't equal each other, so $\|\nabla u^{m,l}\|^2_{L^2(\Omega)}$ $(l=1,\dots,m)$ can't be cancelled out by summation.
I was stuck here and hadn't tried to extend the rest of the original proof to this more general case, so I am not sure whether I will encounter more difficulties as I step further (I think I will because the constants $C$ in my proof by far are dependent on $\varepsilon$ and $h$, which may bring some difficulties when doing approximation). But up to now I hope someone can help me to solve those two questions I raised above, or tell me where I can find any reference. Thank you very much.
