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

## lunedì 20 giugno 2016

## mercoledì 15 giugno 2016

### 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 .

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 .

## mercoledì 1 giugno 2016

### 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 φ

_{1}(A) and PHIO*is**φ*_{0}(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.

Iscriviti a:
Post (Atom)