Finite Difference - Convection Diffusion

This worksheet illustrates the finite difference method for a simple one-dimensional steady convection diffusion equation.

Restart Maple and load libraries.

> restart:

> with(LinearAlgebra):with(plots):

Warning, the name changecoords has been redefined

We consider the equation

-diff(u,x,x)+k*diff(u,x) = 0 u(0) = 0 u(1) = 1

define the parameter k

> k:=50;

k := 50

the number of grid points n and the mesh spacing h

> n:=11;

n := 11

> h:=1/(n-1);

h := 1/10

After discretization we get a system of linear algebraic equations

Au = b

define the matrix and vectors

> A:=Matrix(n,n);

A := _rtable[18266692]

> u:=Vector(n);

u := _rtable[18233348]

> f:=Vector(n);

f := _rtable[16145620]

set f[1] and f[n] to account for boundary conditions

> f[1]:=0;

f[1] := 0

> f[n]:=1;

f[11] := 1

Compute the elements of the matrix A using centered finite differences

> for i from 2 to n-1 do

> A[i,i]:=2/h^2

> end do:

> for i from 2 to n-1 do

> A[i,i-1]:=-k/(2*h)-1/h^2

> end do:

> for i from 2 to n-1 do

> A[i,i+1]:=k/(2*h)-1/h^2

> end do:

set first and last equation to account for boundary condition

> A[1,1]:=1;

A[1,1] := 1

> A[n,n]:=1;

A[11,11] := 1

solve linear system

> u:=LinearSolve(A,f);

u := _rtable[18233428]

Set up for graphics

> q:=[seq([(i-1)*h,u[i]],i=1..n)]:

> ap:=listplot(q,thickness=2,color=black):

Exact solution

> e:=(exp(k*x)-1)/(exp(k)-1);

e := (exp(50*x)-1)/(exp(50)-1)

> ep:=plot(e,x=0..1,thickness=2,linestyle=2,color=black):

Plot and compare approximate solution to exact solution

> display(ep,ap);

[Maple Plot]

Repeat the above (without comments) for a finite difference discretization using first-order

upwind differences for the convective term.

> for i from 2 to n-1 do

> A[i,i]:=2/h^2+k/(2*h)

> end do:

> for i from 2 to n-1 do

> A[i,i-1]:=-k/(2*h)-1/h^2

> end do:

> for i from 2 to n-1 do

> A[i,i+1]:=-1/h^2

> end do:

> A[1,1]:=1;

A[1,1] := 1

> A[n,n]:=1;

A[11,11] := 1

> u:=LinearSolve(A,f);

u := _rtable[18233548]

> q:=[seq([(i-1)*h,u[i]],i=1..n)]:

> ap:=listplot(q,thickness=2,color=black):

> e:=(exp(k*x)-1)/(exp(k)-1);

e := (exp(50*x)-1)/(exp(50)-1)

> ep:=plot(e,x=0..1,thickness=2,linestyle=2,color=black):

> display(ep,ap);

[Maple Plot]

>