Filling points in a grid - Forward Euler algorithm - wrong output

Question

I will very briefly try to explain what I'm doing to those who are less experienced with mathematics, it's really quite simple.

We are trying to fill a grid, as follows:

We find the orange point, U(j,n+1), using three points in a row below it, U(j-1,n), U(j,n), U(j,n+1)

Where the value of U in the entire bottom row is given, and is periodic. So theoretically we can fill this entire grid.

The formula for calculating the orange point is:

U(j,n+1) = U(j,n) + (delta_t / (2 * delta_x)) * (U(j+1,n) - U(j-1,n))

We can write it easily as a system of linear equations as follows:

And now we just repeat this process of multiplying by this matrix (iterating through the time variable) as much as we want. That's a simple way to numerically approximate a solution to a partial differential equation.

I wrote a code that does this, and then I compare my final row, to the known solution of the differential equation.

This is the code

import math
import numpy

def f(x):
    return math.cos(2 * math.pi * x)


def solution(x, t):
    return math.cos(2 * math.pi * (x + t))


# setting everything up
N = 16

Lambda = 10 ** (-20)
Delta_x = 1/(N+1)
Delta_t = Lambda * Delta_x * Delta_x
t_f = 5
v_0 = numpy.zeros((N, 1))

# Filling first row, initial condition was given
for i in range(N):
    v_0[i, 0] = f(i * Delta_x)

# Create coefficient matrix
M = numpy.zeros((N, N))
for i in range(N):
    M[i, i - 1] = -Delta_t / (2 * Delta_x)
    M[i, i] = 1
    M[i, (i + 1) % N] = Delta_t / (2 * Delta_x)

# start iterating through time
v_i = v_0
for i in range(math.floor(t_f / Delta_t) - 1):
    v_i = numpy.dot(M, v_i)

v_final = v_i
if (Delta_t * math.ceil(t_f / Delta_t) != t_f): #we don't reach t_f exactly using Delta_t
    v_final = (1/2) * (v_i + numpy.dot(M, v_i))

u = numpy.zeros(v_final.shape)
for i in range(N):
    u[i, 0] = solution(i * Delta_x, t_f)

for x in range(v_final.shape[0]):
    print (v_final[x], u[x])

theoretically speaking, I should be able to find lambda small enough such that v_final and the known solution, u, will be very similar.

But I can't. No matter how small I make lambda, how finde I make the grid, I seem to converge to something incorrect. They aren't close.

I can't for the life of me figure out the problem. Does anyone have an idea what might be wrong?

Answer 1

You should have Delta_x = 1.0/N, as you divide the interval into N cells.

You get N+1 points on the grid from u[0] to u[N], but as per boundary condition u[N]=u[0], there you also only use an array of length N to hold all the node values.

Per your given formulas you have gamma = dt/(2*dx), thus the reverse computation should be dt = gamma*2*dx or in your variable names

Delta_t = Lambda * 2 * Delta_x

Or you are aiming at the error of the method which is O(dt, dx²) so that it would make sense to have dt = c*dx^2, but not with a ridiculous factor like of c=1e-20, if you want the time discretization error small against the space discretization error, c=0.1 or c=0.01 should be sufficient.

---

import numpy as np

def f(x):
    return np.cos(2 * np.pi * x)

def solution(x, t):
    return f(x + t)


# setting everything up
N_x = 16

Lambda = 1e-2
Delta_x = 1./N_x
Delta_t = Lambda * Delta_x * Delta_x
t_f = 5
N_t = int(t_f/Delta_t+0.5); t_f = N_t*Delta_t

# Filling first row, initial condition was given
x = np.arange(0,N_x,1) * Delta_x
v_0 = f(x)

# Create coefficient matrix
M = np.zeros((N_x, N_x))
for i in range(N_x):
    M[i, i - 1] = -Delta_t / (2 * Delta_x)
    M[i, i] = 1
    M[i, (i + 1) % N_x] = Delta_t / (2 * Delta_x)

# start iterating through time
v_i = v_0[:]
for i in range(N_t):
    v_i = np.dot(M, v_i)

v_final = v_i

u = solution(x, t_f)

for vx, ux in zip(v_final, u):
    print (vx, ux)

---

The Euler method is also not the most precise method, the expected error is in the range exp(L*t_f)*dx^2 = e^5/N_x^2=0.58 for N_x=16 where L=1 was taken as approximate Lipschitz constant. Now if you increase to N_x=50 this error estimate reduces to 0.06 which is also visible in the results.

---

The t exact solution of the x discretized problem is cos(2*pi*(x+c*t)) where c=sin(2*pi*dx)/(2*pi*dx). If you compare against that formula, the errors should be really small of size O(dt).