1
//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?

theSharpShooter
  • 33
  • 8
  • Can you check the signs in the matrix exponentials, currently they look like the original ODE is `u' = A*u + g(t,u)`. Perhaps the result is because you are integrating backward in time. – Lutz Lehmann Mar 31 '17 at 18:51
  • Yes I tried changing the sign also but its answer is still varying for that calculated by Euler method – theSharpShooter Mar 31 '17 at 18:59
  • @LutzL is there something wrong with the implementation? Cuz the formula seems correct? – theSharpShooter Mar 31 '17 at 19:07
  • No, the formula is wrong by a symbol, you switched the factors resp. direction in the division, see answer. – Lutz Lehmann Mar 31 '17 at 19:09
  • https://en.wikipedia.org/wiki/Exponential_integrator I tried to implement equation 2 – theSharpShooter Mar 31 '17 at 19:10
  • Read further down at https://en.wikipedia.org/wiki/Exponential_integrator#First-order_forward_Euler_exponential_integrator. Also note that in Matlab `A/B == A*inv(B)` and `A\B == inv(A)*B`. – Lutz Lehmann Mar 31 '17 at 19:13
  • Your Euler method also has a factor `h` too much, should be `u += h*(A*u + g(t,u))` and corresponds to `u'=Au+g`. – Lutz Lehmann Mar 31 '17 at 19:24
  • What exactly is `Strang`? The tridiagonal matrix of the discretization of the second derivative? – Lutz Lehmann Mar 31 '17 at 19:26
  • Strang(11) gives a matrix in julia that has diagonals = -2 and elements next to diagonals are 1 ie (i+1,i) and (i,i+1) for all except the left and the right corner – theSharpShooter Mar 31 '17 at 19:30
  • Then you are missing the denominator `dx^2` of the differential/difference quotient. The PDE if I read it correctly should be the reaction-diffusion equation `u_t = u_{xx} + (2-u)` which indeed will converge towards the constant `2`. You could also assemble all of `u_{xx}-u` in the linear part. – Lutz Lehmann Mar 31 '17 at 20:09
  • Yes that's a reaction-diffusion and it should converge to 2, But I didn't get the other part, can you please share with me the Matlab code? – theSharpShooter Mar 31 '17 at 22:41

2 Answers2

3

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
Chris Rackauckas
  • 18,645
  • 3
  • 50
  • 81
0

The exponential Euler-Rosenbrock step should be for u'=Lu+g(t,u)

u = expm(h*L)*u + ((expm(h*L)-1)/L)*g(t,u)

Note that the inverted matrix is L (in Matlab A/B == A*inv(B) and A\B == inv(A)*B), and there is no naked h in the formula. In your case, L = -A.

An alternative method to split off the exponential part has is as follows. Set u(t) = exp(-t*A)*v(t) which translates into the differential equation

exp(-t*A)*v'(t) = g(t,exp(-t*A)*v((t))

Now apply the forward Euler step for v

v(t+h) = v(t) + h * exp(t*A)*g(t,exp(-t*A)*v((t))

which translated back into terms in u gives

u(t+h) = exp(-(t+h)*A)*v(t+h) 
       = exp(-h*A) * ( u(t) + h * g(t,u(t)) )

This should equally result in a first order method.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • I was trying to implement exponential rk method, but before that Was thinking of trying exp integrators to solve the same, eq 1.6 of this paper https://na.math.kit.edu/download/papers/acta-final.pdf – theSharpShooter Mar 31 '17 at 19:18
  • Tried your upper part of the answer but its still giving the wrong result!! @LutzL – theSharpShooter Mar 31 '17 at 22:36