Week 4: augmented matrix

In this week I wrote a code that implements the augmented matrix described in theorem 2.1 of [HAM 11]. This matrix will be given as input to expmv and will allow us to implement the exponential integrators by evaluating a single exponential of this augmented matrix avoiding the need to compute any φ functions. The code is

function [Atilde, eta] = augmat(h,A,V)

p = size(V,2);
W = fliplr(V/diag(h.^(0:p-1)));
eta = 2^-ceil(log2(max(norm(W,1),realmin)));
Atilde = [ A , eta*W ; zeros(p,size(A)) , diag(ones(p-1,1),1) ];

in [HAM 11] eta is introduced to avoid rounding errors. I made a small change to avoid further errors, I added the max function so that smaller values ​​than realmin will not be considered.

Now I summarize my work during this month.

Week 1: phi functions. I implemented four functions, based on [BSW 07] (phi1m.m, phi2m.m, phi3m.m, phi4m.m).

Week 2-3: general exponential schemes. I implemented the two schemes for a general exponential Runge-Kutta and Rosenbrock integrator (exprk.m, exprb.m).
As already mentioned, these schemes are not really fast and efficient, but these are working codes that I applied to four different exponential methods (exprk2.m, exprk3.m, exprb3.m, exprb4.m) and I will use them as a reference when I go to implement the official methods.

Week 4: augmented matrix (augmat.m).

My codes can be found here

https://bitbucket.org/ChiaraSegala/chiara-segala

Week 2-3: general exponential schemes

During these two weeks I implemented the two schemes for a general exponential Runge-Kutta and Rosenbrock integrator, see the schemes described in my second post.
I also implemented another file for exponential Runge-Kutta integrators based on the following scheme.
Given the problem
$$
u'(t) = F(t,u) = A u(t) + g(t, u(t)),
\\
u(t_0) = u_0,
$$
the new scheme is
$$
U_{ni} = u_n + c_i h_n \varphi_1(c_i h_n A) F(t_n,u_n) + h_n \sum_{j=2}^{i-1} a_{ij}(h_n A) D_{nj} ,
\\
u_{n+1} = u_n + h_n \varphi_1(h_n A) F(t_n,u_n) + h_n \sum_{i=2}^s b_i(h_n A) D_{ni} ,
$$
where
$$
D_{nj} = g(t_n + c_j h_n, U_{nj}) - g(t_n, u_n).
$$
The main motivation for this reformulation is that the vectors $D_{nj}$ are expected to be small in norm so the application of matrix functions to these vectors are more efficient.

I applied the general pattern to the third-order exponential Runge-Kutta method with tableau

where $\varphi_{j,k} =\varphi_j (c_k h A) $ and $\varphi_j =\varphi_j (h A) $.

Then I also applied the general Rosenbrock pattern to the fourth-order exponential method with a third-order error estimator. Its coefficients are

where $\varphi_j =\varphi_j (h J_n) $.

I tested the correctness and order of the schemes with the following example, a semilinear parabolic problem
$$
\frac{\partial U}{\partial t} (x,t) - \frac{\partial^2 U}{\partial x^2} (x,t) = \frac{1}{1 + U(x,t)^2} + \Phi (x,t)
$$
for $x \in [0,1]$ and $t \in [0,1]$ with homogeneous Dirichlet boundary conditions. $\Phi$ is chosen in such a way that the exact solution of the problem is $U(x,t) = x(1-x)e^{t}$. I discretize this problem in space by standard finite differences with 200 grid points.

For details see [HO 10] and [HO 05].

Within this week I will insert on Bitbucket the codes of phi functions and the three general schemes .

Week 1: phi functions

 

During the first week of GSoC, I wrote four m-files for matrix functions, phi1m,..., phi4m.
I implemented the four functions, based on [BSW 07].
Matrix functions are defined by the recurrence relation

In my files, I use an algorithm based on first a variant of the scaling and squaring approach and after scaling, I use a Padé approximation.

Below the code phi1m.m

function [N, PHI0] = phi1m (A, d)

N is φ₁(A) and PHIO is φ₀(A). d is the order of the Padé approximation.

First, I scale A by power of 2 so that its norm is < ½

l=1;

s= min (ceil (max (0,1+log2 (norm (A, inf)))), 1023);

A = A/2^s;

then I use a (d,d)-Padé approximation

where the polynomials are

I write polynomials in Horner form.

ID = eye (size (A));

i=d;

Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

Dcoeff = ((-1)^(i))*prod ([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

N= Ncoeff;

D= Dcoeff;

for i = (d-1):-1:0

   Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

   N = A*N + Ncoeff * ID ;

   Dcoeff = ((-1)^(i))* prod([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

   D = A*D + Dcoeff * ID ;

endfor

N = full (D\N);

and finally I use the scaling relations

PHI0 = A*N+ID;

for i = 1:s

   N = (PHI0+ID)*N/2;

   PHI0 = PHI0*PHI0;

endfor

in the same way I wrote the other three files.  As soon as possible, I will create a repository on Bitbucket and I will put the codes there.