Scientific Computing I Module 7: Solving the 1D Heat Equation - - PowerPoint PPT Presentation

scientific computing i
SMART_READER_LITE
LIVE PREVIEW

Scientific Computing I Module 7: Solving the 1D Heat Equation - - PowerPoint PPT Presentation

Lehrstuhl Informatik V Scientific Computing I Module 7: Solving the 1D Heat Equation Miriam Mehl based on Slides by Michael Bader (Winter 09/10) Winter 2011/2012 Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing


slide-1
SLIDE 1

Lehrstuhl Informatik V

Scientific Computing I

Module 7: Solving the 1D Heat Equation

Miriam Mehl based on Slides by Michael Bader (Winter 09/10)

Winter 2011/2012

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 1

slide-2
SLIDE 2

Lehrstuhl Informatik V

Part I: Analytical Solution

The Heat Equation in 1D Analytical Solutions Computing Analytical Solutions – Steps Separation of Variables Ensuring the Initial Conditions – Fourier’s Method

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 2

slide-3
SLIDE 3

Lehrstuhl Informatik V

The Heat Equation in 1D

  • remember the heat equation without heat sources:

Tt = κ∆T

  • we examine the 1D case, and set κ = 1 to get:

ut = uxx for x ∈ (0, 1), t > 0

  • using the following initial and boundary conditions:

u(x, 0) = g(x), x ∈ (0, 1) u(0, t) = u(1, t) = 0, t > 0

t x

?

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 3

slide-4
SLIDE 4

Lehrstuhl Informatik V

Computing Analytical Solutions – Steps

Steps:

  • 1. try to find some solution of the PDE
  • 2. try to satisfy boundary conditions
  • 3. try to satisfy initial conditions

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 4

slide-5
SLIDE 5

Lehrstuhl Informatik V

Separation of Variables

  • assumption:

u(x, t) = X(x) · T(t)

  • insert this assumption into the heat equation

∂ ∂t (X(x) · T(t)) = ∂2 ∂x2 (X(x) · T(t))

  • r

X(x) · Tt(t) = T(t) · Xxx(x)

  • divide by X(x)T(t), and get:

Tt(t) T(t) = Xxx(x) X(x)

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 5

slide-6
SLIDE 6

Lehrstuhl Informatik V

Transforming the PDE into two ODEs

  • separation of variables leads to:

Tt(t) T(t) = Xxx(x) X(x) = −λ.

  • thus, we obtain two ODEs:

Xxx(x) + λX(x) = (1) Tt(t) + λT(t) = 0. (2)

  • solve X(x)-part:

X(x) = sin √ λx

  • r X(x) = cos

√ λx

  • ,
  • solve T(t)-part:

T(t) = e−λt.

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 6

slide-7
SLIDE 7

Lehrstuhl Informatik V

Ensuring the Boundary Conditions

  • We have computed general solutions

u(x, t) = sin √ λx

  • e−λt or u(x, t) = cos

√ λx

  • e−λt.

0.02 0.04 0.06 0.08 0.1 0.5 1 −1 −0.5 0.5 1 t x

The boundary conditions of

  • ur example are

u(0, t) = u(1, t) = 0 for all t.

  • Thus, we can reduce possible solutions to

uk(x, t) = sin (kπx) e−(kπ)2t, k = 1, 2, 3 . . .

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 7

slide-8
SLIDE 8

Lehrstuhl Informatik V

Ensuring the Initial Conditions – Fourier’s Method

  • The functions

uk(x, t) = sin(kπx)e−(kπ)2t, k = 1, 2, . . . , solve the 1D heat equation PDE with correct boundary values and for the initial conditions: uk(x, 0) = sin(kπx), x ∈ (0, 1).

0.02 0.04 0.06 0.08 0.1 0.5 1 0.2 0.4 0.6 0.8 1 t x 0.02 0.04 0.06 0.08 0.1 0.5 1 −1 −0.5 0.5 1 t x 0.02 0.04 0.06 0.08 0.1 0.5 1 −1 −0.5 0.5 1 t x

  • To ensure the initial conditions, we use the Fourier sine series:

g(x) =

  • k=1

ck sin(kπx)

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 8

slide-9
SLIDE 9

Lehrstuhl Informatik V

Limits of Analytical Solutions

  • For the 1D heat equation, we obtain an analytical solution

u(x, t) =

  • k=1

cke−(kπ)2t sin(kπx).

  • Solutions for the heat equation with heat sources?

Separation of variables yields X(x) · Tt(t) = T(t) · Xxx(x) + f(x). ⇒ No analytical solution computable for arbitrary f.

  • need for numerical methods!!!

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 9

slide-10
SLIDE 10

Lehrstuhl Informatik V

Part II: Numerical Solution

Numerical Solution 1 – An Explicit Scheme Discretisation Accuracy of the Explicit Scheme Stability of the Explicit Scheme Numerical Solution 2 – An Implicit Scheme Implicit Time-Stepping Stability of the Implicit Scheme

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 10

slide-11
SLIDE 11

Lehrstuhl Informatik V

Numerical Solution 1 – An Explicit Scheme

Discretisation similar to ODEs:

  • compute approximations

v(m)

j

≈ u(xj, tm) at grid points xj and time points tm: xj := j · h tm := m · τ,

t

t1 t2 x

1

x

2

x

N

tM

x

? approximate uxx at (xj, tm) by v(m)

j−1 − 2v(m) j

+ v(m)

j+1

h2 , ut at (xj, tm) by v(m+1)

j

− v(m)

j

τ .

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 11

slide-12
SLIDE 12

Lehrstuhl Informatik V

An Explicit Scheme

  • This approximation results in the equations

v(m+1)

j

− v(m)

j

τ = v(m)

j−1 − 2v(m) j

+ v(m)

j+1

h2 . (3) for j = 1, . . . , N − 1 and m ≥ 0.

  • Add initial and boundary conditions:

v(m) = v(m)

n

= 0, for all m ≥ 0, v(0)

j

= g(xj), for j = 1, . . . , n − 1

  • and obtain an explicit scheme:

v(m+1)

j

= v(m)

j

+ τ h2

  • v(m)

j−1 − 2v(m) j

+ v(m)

j+1

  • .
  • We can, step by step, compute v(m)

j

starting with v(0)

j

= g(xj).

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 12

slide-13
SLIDE 13

Lehrstuhl Informatik V

Accuracy of the Explicit Scheme

Observations:

  • first order accurate in τ:

ut(x, t) = u(x, t + τ) − u(x, t) τ + O(τ)

  • second order accurate in h:

uxx(x, t) = u(x − h, t) − 2u(x, t) + u(x + h, t) h2 + O(h2)

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 13

slide-14
SLIDE 14

Lehrstuhl Informatik V

Stability of the Explicit Scheme

  • exact solution:

uk(x, t) = e−(kπ)2t sin(kπx)

  • explicit scheme:

v(m+1)

j

− v(m)

j

τ = v(m)

j−1 − 2v(m) j

+ v(m)

j+1

h2 Is there a stability condition for τ as observed in the ODE case?

  • assumption: there are solutions of the form

v(m)

j

= (ak)m sin(πkxj), where xj := jh. (4)

  • compare with the exact solution ⇒ (ak)m should decrease

⇒ |ak| < 1 necessary

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 14

slide-15
SLIDE 15

Lehrstuhl Informatik V

Stability of the Explicit Scheme (2)

  • Inserting v(m)

j

= (ak)m sin(πkxj) into v(m+1)

j

− v(m)

j

τ = v(m)

j−1 − 2v(m) j

+ v(m)

j+1

h2 leads to ak = 1 + τ h2 (2 cos(πkh) − 2) = 1 − 4τ h2 sin2 πkh 2

  • .
  • for stability: |ak| < 1 for all k if

4τ h2 < 2

  • r

τ < h2 2 .

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 15

slide-16
SLIDE 16

Lehrstuhl Informatik V

Stability of the Explicit Scheme (3)

  • solution for general initial conditions:

v(0)

j

= fj =

n−1

  • k=1

ck sin(πkxj), v(m)

j

:=

n−1

  • k=1

ck(ak)m sin(πkxj)

  • stability for τ < h2

2 because |ak| < 1 for all k

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 16

slide-17
SLIDE 17

Lehrstuhl Informatik V

Numerical Solution 2 – An Implicit Scheme

  • apply implicit Euler:

v(m+1)

j

− v(m)

j

τ = v(m+1)

j−1

− 2v(m+1)

j

+ v(m+1)

j+1

h2 for j = 1, . . . , n − 1, and m ≥ 0.

  • boundary conditions:

v(m) = v(m)

n

= 0, for all m ≥ 0,

  • initial conditions:

v(0)

j

= g(xj), for j = 1, . . . , n − 1.

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 17

slide-18
SLIDE 18

Lehrstuhl Informatik V

An Implicit Scheme

  • implicit equations for v(m+1)

j

: v(m+1)

j

− τ h2

  • v(m+1)

j−1

− 2v(m+1)

j

+ v(m+1)

j+1

  • = v(m)

j

.

  • With the ratio r := τ/h2, we can write it as

−rv(m+1)

j−1

+ (1 + 2r)v(m+1)

j

− rv(m+1)

j+1

= v(m)

j

for j = 1, . . . , n − 1, and m ≥ 0.

  • solution:

v(m) = (I + rA)−1v(m−1), where A = tridiag(−1, 2, −1).

  • Use a solver for systems of linear equations to obtain v(m+1)

j

in every step.

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 18

slide-19
SLIDE 19

Lehrstuhl Informatik V

Stability of the Implicit Scheme

  • implicit scheme:

v(m+1)

j

− v(m)

j

τ = v(m+1)

j−1

− 2v(m+1)

j

+ v(m+1)

j+1

h2

  • again: insert assumed solutions

u(m)

j

= (ak)m sin(πkxj), where xj := jh

  • into the implicit scheme, and obtain:

ak =

  • 1 + 4τ

h2 sin2 πkh 2 −1 .

  • 0 < ak < 1 independent of k and h

⇒ unconditionally stable

Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 7: Solving the 1D Heat Equation, Winter 2011/2012 19