So here's the background on the problem at large, it won't be useful I think for what I'm asking help on since the formulation has nothing to do with finding eigenvalues. But it might be useful in giving clarity to the situation.

So consider applying the method of lines to the heat equation


with initial conditions $u(0,x)=g(x),\, 0\leq x\leq 1$ and boundary conditions $u(t,0)=u(t,1)=0,\, t\geq 0$.

For simplicity, assume $\Delta x=\frac{1}{m+1}$ and let $y_i(t)$ approximate $u(x_i,t)$ where $x_i=i\Delta x,\, i=0,1,...,m+1$.

Then using second-order difference scheme to approximate the heat equation we arrive at

$\frac{dy_i}{dt}=\frac{a(y_{i+1}-2y_i+y_{i-1})}{(\Delta x)^2},\, i=1,...,m$ with $y_0(t)=0$ and $y_{m+1}(t)=0$.

So now we get a system of ODE's with initial data $\vec{y}'=A\vec{y}$ where $A=\frac{a}{(\Delta x)^2}\begin{bmatrix} -2&1&0&.&.&0 \\ 1&-2&1&.&.&.\\ 0&1&-2&1&.&.\\.&.&.&.&.&. \\.&.&.&1&-2&1\\0&.&.&.&1&-2 \end{bmatrix}$

(Sorry if my matrix looks half-assed, but this is supposed to represent a symmetric block diagonal matrix)

FYI: The component form for $\vec{y}'$ and $\vec{y}$ are trivial from numerical scheme.

Now, for this next part we do a bit hand wavy since the focus isn't so much solving ODE's or PDE's the natural way like using separation of variables; we are focusing more on the numerical methods for them, so we suppose we are given the eigenfunctions beforehand, and they are,

$\vec{v}^k=(v_1^k,v_2^k,...,v_m^k)^T=(\sin(k\pi\Delta x),\sin(2k\pi\Delta x),...,\sin(mk\pi\Delta x))^T,\, 1\leq k\leq m$, or in simplified form $v_i^k=\sin(ik\pi\Delta x),\, i=1,...,m$

Now to my first question, how do I derive the eigenvalues $\lambda_k$ without making a huge mess making busy calculation? cause if it does involve that then I'd get on it and won't need your help on this part, but providing a more simple way of doing this that involves ideas discussed in numerical schemes for these equations is what I should be looking to do. But lets say we carried out this brute force calculation of the eigenvalues the natural way, wouldn't this lead to different eigenvalues for each row since they all have different trig functions?

I probably could've asked this question without much of the PDE and numerical method jargon used to derive the matrix, but I wanted to show that I actually know what I'm trying to do.

If you can provide me with the eigenvalue that would be all that is needed, I can figure out how to get there from what I have given you so far, unless it is super trivial and is a one-line calculation.

Now, to my second question, how would I use the eigenvalues of A to show that the forward Euler discretization of the method of lines is stable if $h<\frac{1}{2a}(\Delta x)^2$ Why would we need a stability restriction for the time step for the forward Euler but not the backward Euler(Crank-Nicolson method)?

Someone asked this similar problem but solutions were insufficient and didn't provide enough clarity of the concepts being discussed

Method of Lines Diffusion Problem

  • $\begingroup$ Since I posted this, I've come to the conclusion that determining the eigenvalues might be solvable using a method for ODE systems with constant coefficients. Since A is constant, diagonalizable, $m\times m$ matrix then there exists a nonsingular matrix T, consisting of the eigenvectors of A such that $T^{-1}AT=\Lambda$. Where $\Lambda=diag\{\lambda_1,\lambda_2,...,\lambda_m\}$. Now it makes sense as to why the problem gave us the form $v_i^k=\sin(ik\pi\Delta x),\, i=1,...,m,\, 1\leq k\leq m$ $\endgroup$ Nov 1 '21 at 15:56

For an eigen-solution one expects that $y_i(t)=c(t)v_i$ where here $v_i=\sin(\omega x_i)$, and $ω$ getting its discrete values from $v_{M+1}=0$. This you can insert directly in the discretization equation without going to the matrix formulation $$ c'(t)v_i=\frac{a}{Δx^2}c(t)(v_{i+1}-2v_i+v_{i-1}). $$ The last term can now be attacked using trigonometric identities, like $$ \sin(A+B)+\sin(A-B)=2\sin(A)\cos(B), $$ to get $$ c'(t)v_i=\frac{a}{Δx^2}c(t)(2\cos(ωΔx)-2)v_i. $$ Thus the relevant eigenvalue is $\lambda=-\frac{4a\sin^2(ωΔx/2)}{Δx^2}$. As the intersection of the stability region and the real line for explicit Euler is the interval $[-2,0]$, stability is satisfied for $$ 2\ge\frac{4ah}{Δx^2}\ge λh. $$ For backwards Euler the stability region contains $(-\infty,0]$, so that no restrictions on the step size result from this side.

  • $\begingroup$ where are you getting $c(t)$ from? Is there an explicit form that it takes from all the variables given? $\endgroup$ Nov 1 '21 at 15:47
  • $\begingroup$ I have added an addendum to my post under comments, will this work as well in giving the eigenvalue you obtained? $\endgroup$ Nov 1 '21 at 16:04
  • $\begingroup$ You transformed $u_t=au_{xx}$ to $y'=Ay$. Now if $v$ is an eigenvector of $A$, $Av=λv$, any solution that starts in direction $v$ remains in that direction, so $y(t)=c(t)v$ where $c(t)$ scalar with $c'(t)=λc(t)$, which is now easily solvable as exponential function. Also the action of the Euler or other Runge-Kutta step is similarly simplified in the direction of eigenvectors. $\endgroup$ Nov 1 '21 at 16:24
  • $\begingroup$ Oh, I see. It's like that test function that you use to show whether or not your numerical method is stable or not, right? Depending on where in the Complex plane you solution falls in $\endgroup$ Nov 1 '21 at 16:34
  • $\begingroup$ Yes, in the eigenbasis the system decouples into a system of scalar "test" equations for the different eigenvalues. $\endgroup$ Nov 1 '21 at 16:44

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.