Why is exponential integrator method giving poor results?

Question

//u' + Au = g(t,u) can be solved by exponential integrators also
//Following snippet is for exp INtegrators
A = -full(Strang(11)) 
A[end,1]=1;A[1,end]=1;
g(t,u) = 2-u
u0 = zeros(11);u0[6]=1
nsteps = 1000
tmax = 10.0
h = tmax/nsteps
u = u0
t = 0
for k in 1:nsteps
    u = expm(-h*A)*u + h*((expm(-h*A)-1)\(-h*A))*g(t,u)
    t = k*h
end
//this is for euler's method
for k in 1:nsteps
    u += h*(A*u + h*g(t,u))
    t = k*h
end

Why are they giving poor results?

The method is exploding very badly, it should converge to [1.99]*11 , or something like that?

Is there any mistake while implementing Exp Integrator?

Answer 1

The test problem is a singular matrix. A better test is the setup:

using SpecialMatrices
A = -full(Strang(11))
g(t,u) = 2-u
u = zeros(11);u[6]=1
nsteps = 10000
tmax = 1.0
h = tmax/nsteps
t = 0

Using this, fix the h in the Euler to get (notice there's an extra h, my bad:

u = zeros(11);u[6]=1
for k in 1:nsteps
    u += h*(A*u + g(t,u))
    t = k*h
end
@show u

u = [0.93573,1.19361,1.26091,1.29627,1.34313,1.37767,1.34313,1.29627,1.26091,1.19361,0.93573]

But to find out what's wrong, start looking at numbers. What happens A=0? Well, we know that phi(z) = (e^z - 1)/z. By L'Hopital's rule, phi(z) -> 1 as z->0. Therefore, in order for our implementation to have the same behavior, we have to have that same result. But let's check what happens:

expm(zeros(5,5))

5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

Notice that this gives the identity matrix. So think about the limit: if the bottom is going to zero... how can this be constant? We must have that the top is going to zero... so the top is going to I.

And that's the moment of clarity: the author meant 1 in the field that you're in. So for a matrix input, 1=I. When you realize that, you fix the code:

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = expm(h*A)*u + ((expm(h*A)-I)/A)*g(t,u)
  t = k*h
end
@show u

u = [0.935722,1.1936,1.26091,1.29627,1.34313,1.37768,1.34313,1.29627,1.26091,1.1936,0.935722]

Moral of the story: for programming mathematics, you also have to debug your math.

Edit

Get a more efficient form one step at a time. First, try and force another varphi term:

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = (I + A*(expm(h*A)-I)/A)*u + ((expm(h*A)-I)/A)*g(t,u)
  t = k*h
end
@show u

Now gather:

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = u + ((expm(h*A)-I)/A)*(A*u + g(t,u))
  t = k*h
end
@show u

This is the efficient form of the method you were trying to write. Then you can cache the operator since A is constant:

# Norsett-Euler
u = zeros(11);u[6]=1
phi1 = ((expm(h*A)-I)/A)
for k in 1:nsteps
  u = u + phi1*(A*u + g(t,u))
  t = k*h
end
@show u