Heat and Diffusion Equation Schemes

( Heat Equation and Stability of the Scheme). Consider the implicit in time, Backward Euler method for the solution of the heat equation: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} -10\frac {u_{m+1}^{n+1}-2u_m^{n+1}+u_{m-1}^{n+1}}{h^2} = f_m^{n+1}, \] \[ u_m^0 = g_m,\qquad m = 0,\pm 1,\pm 2,\ldots ,\qquad n = 0,1,\ldots ,[T/\Delta t]-1, \]

and investigate the stability of the scheme using the von Neumann analysis.

Heat Equation Stability

Consider the variable coefficient diffusion equation \[ v_t=(\beta (x)v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

DO NOT NEGLECT THE FACT THAT THE PROBLEM HAS VARIABLE COEFFICIENTS AND THAT THERE ARE BOUNDARY CONDITIONS AT 0 AND 1!

Heat Equation Stability:

Consider the variable coefficient diffusion equation \[ v_t=(\beta (x)v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

DO NOT NEGLECT THE FACT THAT THE PROBLEM HAS VARIABLE COEFFICIENTS AND THAT THERE ARE BOUNDARY CONDITIONS AT 0 AND 1!

Fact 1: A real symmetric \(n\times n\) matrix \(A\) can be diagonalized by an orthogonal similarity transformation, and \(A\)’s eigenvalues are real.

Fact 2: The \((N-1)\times (N-1)\) matrix \(M\) defined by \[ M= \begin {pmatrix} -2 & 1 & 0 & \cdots & 0\\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & -2 & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 1\\ 0 & \cdots & 0 & 1 & -2 \end {pmatrix} \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2N}\right ), \qquad l=1,2,\ldots ,N-1. \]

Fact 3: The \((N+1)\times (N+1)\) matrix: \[ \begin {pmatrix} -1 & 1 & 0 & \cdots & 0\\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & -2 & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 1\\ 0 & \cdots & 0 & 1 & -1 \end {pmatrix} \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2(N+1)}\right ), \qquad l=0,1,\ldots ,N. \]

Fact 4: For a real \(n\times n\) matrix \(A\), the Rayleigh quotient of a vector \(x\in \mathbb {R}^n\) is the scalar \[ r(x)=\frac {x^T Ax}{x^T x}. \] The gradient of \(r(x)\) is \[ \nabla r(x)=\frac {2}{x^T x}(Ax-r(x)x). \] If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and \(\nabla r(x)=0\).

Heat Equation Stability:

Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\geq \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

NOTE: Because of the boundary conditions and the variable coefficients, von Neumann (Fourier) analysis does not apply.

7) Numerical Methods for ODEs: Consider the initial value problem for a system of ODEs \[ y'=Ay \] where \(A\) has eigenvalues with real parts ranging from \(-10^6\) to \(-1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. Does your answer change depending on the size of the system and the condition number of \(A\)? Explain your answer. \[ y^{n+1}=y^n+kAy^n, \tag {1} \] \[ y^{n+1}=y^n+kAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {k}{2}\left (Ay^n+Ay^{n+1}\right ). \tag {3} \]

Heat Equation Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0.0,\qquad v(1,t)=0.0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. For any fixed value of \(\theta \in [0,1]\), consider the \(\theta \)-scheme for this problem: \[ \begin {aligned} \frac {u_j^{n+1}-u_j^n}{k} =\frac {1}{h^2}\{& (1-\theta )(\beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1})\\ &+\theta (\beta _{j-1/2}u_{j-1}^n -(\beta _{j-1/2}+\beta _{j+1/2})u_j^n +\beta _{j+1/2}u_{j+1}^n)\}. \end {aligned} \] with appropriate treatment of the boundary conditions in the relevant equations. Here, \(h\) is the spacestep, \((J+1)h=1\), \(k\) is the timestep, and \(\beta _{j+1/2}=\beta (x_{j+1/2})\), where \(x_{j+1/2}=(j+1/2)h\). Let \(M\) be the symmetric \(J\times J\) matrix defined by \[ \setlength {\arraycolsep }{3pt} M= \begin {pmatrix} -(\beta _{1/2}+\beta _{3/2}) & \beta _{3/2} & \cdots & \cdots & \cdots & \cdots \\ \beta _{3/2} & -(\beta _{3/2}+\beta _{5/2}) & \beta _{5/2} & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \beta _{J-3/2} & -(\beta _{J-3/2}+\beta _{J-1/2}) & \beta _{J-1/2}\\ \cdots & \cdots & \cdots & \cdots & \beta _{J-1/2} & -(\beta _{J-1/2}+\beta _{J+1/2}) \end {pmatrix}. \] Analyze the stability of this scheme in terms of \(\theta \), \(\frac {k}{h^2}\), and the eigenvalues of \(M\). Are there values of \(\theta \) for which the scheme is stable for all choices of \(\frac {k}{h^2}>0\)? (Note that answering this question should not require lots of writing.)

1. Derive the “diffusion” or “heat” equation in several spatial dimensions.
2. In the derivation of this equation, it is routinely assumed that \[ \frac {d}{dt}\int _\Omega u(x,t)\,dx=\int _\Omega \frac {\partial u}{\partial t}\,dx. \] Under what conditions is this interchange of order valid? Prove your result. Hint: Examine the function \[ w(t)=\int _0^t\left (\int _\Omega \frac {\partial u(x,t)}{\partial t}\,dx\right )\,dt. \]

Describe quantitatively the diffusion of 3 units of a diffusing material which is placed at the initial moment \(t=0\) at the origin of an infinite two-dimensional plane. Assume that the plane has unit conductivity.