1
$\begingroup$

I am solving the heat equation discretized spatially via FEM and temporally via backward Euler. I get the system $$M \dot{u} = K u +f$$ where $u$ is a vector representing the solution at spatial locations.

For piecewise linear basis functions, the mass matrix is defined as enter image description here while the stiffness matrix $K$ is enter image description here

By discretizing this in time, I get $$u_{k+1} = Au_k + By_k$$ for some matrix $B.$ Here, $A = (I-\Delta t M^{-1}K)^{-1}$ where $\Delta t$ is the time step size.

Is it possible to show that $$\|A\|_2 < 1$$ for $\Delta t$ sufficiently small? I would like an analytical solution without computing this in Matlab.

I have a hard time analyzing this because $M^{-1}K$ is not symmetric although $M,K$ are symmetric. Otherwise, I could have just used the fact that the spectral radius and the 2 norm coincide.

Improve this question
$\endgroup$
1
$\begingroup$

For stability analysis, we can assume $f = 0$. First we write $$ M \left( u_{k+1} - u_{k} \right) = \Delta t K u_{k+1} $$ or $$ \left( M - \Delta t K \right) u_{k+1} = M u_{k} $$ Take the eigendecomposition of $(K, M)$ (there is an explicit expression in 1D), such that $$ w_i^T K w_j = \delta_{ij} \lambda_i \quad \mbox{and} \quad w_i^T M w_j = \delta_{ij} $$ Write $$ w_i^T \left( M - \Delta t K \right) u_{k+1} = w_i^T M u_{k} $$ and express $u_k$ and $u_{k+1}$ in the basis of eigenmodes. From there, you should be able to get an analytical expression to control $\Delta t$.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Hi, thanks for the reply. I am confused about 2 things: 1) I'm not interested in showing stability but rather I want to show that $\|A\|_2 < 1$. Does the approach you mentioned above work for this or is it only to show that the spectral radius is less than 1? 2) What do you mean by express $u_k,u_{k+1}$ as eigenmodes? I have never heard of the method you showed above before. Do you have a reference for this? $\endgroup$ – user1237300 May 1 '20 at 1:49
  • $\begingroup$ 1. Why do you want to show that $||A||_2 < 1$? --- 2. $u_k = \sum_{j = 1}^{n} \mu_k^j w_j$. Using the orthogonality among the eigenvectors $(w_i)$, we get an expression between $\mu_{k+1}^{j}$ and $\mu_{k}^{j}$. $\endgroup$ – user7440 May 1 '20 at 1:56
  • $\begingroup$ Note also $|| ( I - M^{-1} K )^{-1} x ||_2 = || ( M^{-1/2} ( I - M^{-1/2} K M^{-1/2} ) M^{+1/2} )^{-1} x ||_2 = || M^{-1/2} ( I - M^{-1/2} K M^{-1/2} )^{-1} M^{+1/2} x ||_2 = || ( I - M^{-1/2} K M^{-1/2} )^{-1} M^{+1/2} x ||_{M^{-1}}$ $\endgroup$ – user7440 May 1 '20 at 2:00
  • $\begingroup$ I'm using it for error estimates. How does the above help in finding the 2-norm of $I-\Delta t M^{-1}K$? $\endgroup$ – user1237300 May 1 '20 at 2:08
  • $\begingroup$ When you use the expression with $M^{-1}$-norm, the operator $I - \Delta t M^{-1/2} K M^{-1/2}$ is now symmetric and its eigenvectors are the ones of $(K, M)$. So studying $|| ( I - \Delta t M^{-1} K )^{-1} x ||_2 / || x ||_{2}$ is the same as studying $|| ( I - \Delta t M^{-1/2} K M^{-1/2})^{-1} y ||_{M^{-1}} / || y ||_{M^{-1}}$ (where $y = M^{1/2} x$) $\endgroup$ – user7440 May 1 '20 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.