Euler and Backward Euler Time Discretizations

This worksheet illustrates the Euler and Backward Euler methods for time discretization for a simple one-dimensional unsteady convection diffusion equation.

Restart Maple and load libraries.

> restart:

> with(LinearAlgebra):with(plots):

We consider the equation

diff(u,t)-1/5*diff(u,x,x)+4/5*diff(u,x) = 0 u(0,t) = 0 u(Pi,t) = 0

with initial condition

u(x,0) = exp(2*x)*sin(x)

the exact solution for this equation is

u(x,t) = exp(2*x-t)*sin(x)

Set the parameters... things (notation) should be similar to previous Maple worksheet

> k:=4/5;

k := 4/5

> nu:=1/5;

nu := 1/5

the number of grid points n and the mesh spacing h

> n:=41;

n := 41

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

h := 1/40*Pi

the time step k

> deltaT:=0.05;

deltaT := .5e-1

construct the matrix and vectors (unknowns and r.h.s.) as in previous example

> A:=Matrix(n,n);

A := _rtable[15978364]

> f:=Vector(n);

f := _rtable[18232892]

> f[1]:=0;

f[1] := 0

> f[n]:=0;

f[41] := 0

> for i from 2 to n-1 do

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

> end do:

> for i from 2 to n-1 do

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

> end do:

> for i from 2 to n-1 do

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

> end do:

> Id:=IdentityMatrix(n,n);

Id := _rtable[16175844]

The matrix A above is the matrix obtained from semi-discretization, and the matrix Ap below is the one for the fully discrete system

> Ap:=MatrixAdd(Id,A,1,-deltaT);

Ap := _rtable[18276824]

set initial conditions

> u0:=Vector(n);

u0 := _rtable[18233092]

> for i from 2 to n-1 do

> u0[i]:=exp(2*(i-1)*h)*sin((i-1)*h)

> end do:

The following are 5 time steps for the explicit method (Forward Euler)

> u1:=MatrixVectorMultiply(Ap,u0);

u1 := _rtable[18233172]

> u2:=MatrixVectorMultiply(Ap,u1);

u2 := _rtable[18233292]

> u3:=MatrixVectorMultiply(Ap,u2);

u3 := _rtable[18233372]

> u4:=MatrixVectorMultiply(Ap,u3);

u4 := _rtable[18233492]

> u5:=MatrixVectorMultiply(Ap,u4);

u5 := _rtable[18233612]

set up vectors for graphics to display the exact and approximate solutions

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

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

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

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

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

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

> pq0:=listplot(q0,thickness=2,color=black):

> ex0:=plot(exp(2*x)*sin(x),x=0..Pi):

> pq1:=listplot(q1,thickness=2,color=black):

> ex1:=plot(exp(2*x-deltaT)*sin(x),x=0..Pi):

> pq2:=listplot(q2,thickness=2,color=black):

> ex2:=plot(exp(2*x-2*deltaT)*sin(x),x=0..Pi):

> pq3:=listplot(q3,thickness=2,color=black):

> ex3:=plot(exp(2*x-3*deltaT)*sin(x),x=0..Pi):

> pq4:=listplot(q4,thickness=2,color=black):

> ex4:=plot(exp(2*x-4*deltaT)*sin(x),x=0..Pi):

> pq5:=listplot(q5,thickness=2,color=black):

> ex5:=plot(exp(2*x-5*deltaT)*sin(x),x=0..Pi):

> display(pq0,ex0);

[Maple Plot]

> display(pq1,ex1);

[Maple Plot]

> display(pq2,ex2);

[Maple Plot]

> display(pq3,ex3);

[Maple Plot]

> display(pq4,ex4);

[Maple Plot]

> display(pq5,ex5);

[Maple Plot]

Now do similar computations for the implicit method (Backward Euler)

> deltaT:=0.05;

deltaT := .5e-1

> Ap:=MatrixAdd(Id,A,1,+deltaT);

Ap := _rtable[17052684]

> u0:=Vector(n);

u0 := _rtable[18233812]

> for i from 2 to n-1 do

> u0[i]:=exp(2*(i-1)*h)*sin((i-1)*h)

> end do:

> u1:=LinearSolve(Ap,u0);

u1 := _rtable[18233932]

> u2:=LinearSolve(Ap,u1);

u2 := _rtable[18234012]

> u3:=LinearSolve(Ap,u2);

u3 := _rtable[18234092]

> u4:=LinearSolve(Ap,u3);

u4 := _rtable[18234172]

> u5:=LinearSolve(Ap,u4);

u5 := _rtable[18234252]

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

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

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

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

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

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

> pq0:=listplot(q0,thickness=2,color=black):

> ex0:=plot(exp(2*x)*sin(x),x=0..Pi):

> pq1:=listplot(q1,thickness=2,color=black):

> ex1:=plot(exp(2*x-deltaT)*sin(x),x=0..Pi):

> pq2:=listplot(q2,thickness=2,color=black):

> ex2:=plot(exp(2*x-2*deltaT)*sin(x),x=0..Pi):

> pq3:=listplot(q3,thickness=2,color=black):

> ex3:=plot(exp(2*x-3*deltaT)*sin(x),x=0..Pi):

> pq4:=listplot(q4,thickness=2,color=black):

> ex4:=plot(exp(2*x-4*deltaT)*sin(x),x=0..Pi):

> pq5:=listplot(q5,thickness=2,color=black):

> ex5:=plot(exp(2*x-5*deltaT)*sin(x),x=0..Pi):

> display(pq0,ex0);

[Maple Plot]

> display(pq1,ex1);

[Maple Plot]

> display(pq2,ex2);

[Maple Plot]

> display(pq3,ex3);

[Maple Plot]

> display(pq4,ex4);

[Maple Plot]

> display(pq5,ex5);

[Maple Plot]