Difference methods for Parabolic Equations Algorithms

1. Difference methods for Parabolic Equations Algorithms

2. Explicit Dirichlet Preliminaries: define variables: real a2, diffusivity real L, length of domain real k, time step (Dt) real h, space step (Dx) integer jmax, number of time steps integer nmax, number of grid points define functions (may call subroutines) S(x,t), source term f(x), initial condition p(t), Dirichlet boundary condition at x=0 q(t), Dirichlet boundary condition at x=L read input data (a2, L, …)

3. Explicit Dirichlet Define grid: h = L / (nmax + 1) r = a2k/h2 Initialization V(0) = p(0) = f(0) (these should be consistent) V(nmax) = q(0) = f(nmax+1) (these should be consistent) for n = 1, 2, …, nmax V(n) = f(n) endloop

4. Explicit Dirichlet Time loop for j = 1, 2, …, jmax for n = 1, 2, …, nmax U(n) = r*V(n-1)+(1-2*r)*V(n)+r*V(n+1)+k*S(n,j) endloop U(0) = p(j) U(nmax+1) = q(j) Output U Prepare for next time step for n = 0, 1, …, nmax+1 V(n) = U(n) endloop endloop

5. Tridiagonal Matrix

6. Thomas Algorithm Preliminaries define variables Decomposition for n = 2, 3, …, N a(n) = a(n)/b(n-1) b(n) = b(n) – a(n)*c(n-1) endloop Forward Substitution for n = 2, 3, …, N d(n) = d(n) – a(n)*d(n-1) endloop Backward Substitution u(N) = d(N)/b(N) for n = N-1, N-2, …, 1 u(n) = (d(n) – c(n)*u(n+1))/b(n) endloop

7. Implicit Algorithm Preliminaries: define variables: real a2, diffusivity real L, length of domain real k, time step (Dt) real h, space step (Dx) integer jmax, number of time steps integer nmax, number of grid points define functions (may call subroutines) S(x,t), source term f(x), initial condition p(t), Dirichlet boundary condition at x=0 q(t), Dirichlet boundary condition at x=L read input data (a2, L, …)

8. Implicit Algorithm Define grid: h = L / (nmax + 1) r = a2k/h2 Initialization U(0) = p(0) = f(0) (these should be consistent) U(nmax) = q(0) = f(nmax+1) (these should be consistent) for n = 1, 2, …, nmax U(n) = f(n) endloop

9. Implicit Algorithm • Time loop • for j = 1, 2, …, jmax • for n = 1, nmax • a(n) = -r • b(n) = 1+2*r • c(n) = -r • d(n) = U(n) + k*S(n,j) • endloop • d(1) = d(1) + r*p(j) • d(nmax) = d(nmax) + r*q(j) • call tridiagonal solver (nmax, a, b, c, d, U) • U(0) = p(j) • U(nmax+1) = q(j) • endloop

10. Neumann BCs

11. Neumann BCs Preliminaries: define variables: real a2, diffusivity real L, length of domain real k, time step (Dt) real h, space step (Dx) integer jmax, number of time steps integer nmax, number of grid points define functions (may call subroutines) S(x,t), source term f(x), initial condition p(t), Neumann boundary condition at x=0 q(t), Neumann boundary condition at x=L read input data (a2, L, …)

12. Neumann BCs Define grid: h = L / (nmax - 1) r = a2k/h2 Initialization for n = 1, 2, …, nmax U(n) = f(n) endloop

13. Neumann BCs • Time loop • for j = 1, 2, …, jmax • for n = 1, 2, …, nmax • a(n) = -r • b(n) = 1+2*r • c(n) = -r • d(n) = U(n) + k*S(n,j) • endloop • d(1) = d(1) – 2*r*h*p(j) • d(nmax) = d(nmax) + 2*r*h*q(j) • c(1) = -2*r • a(nmax) = -2*r • call tridiagonal solver (nmax, a, b, c, d, U) • endloop