# Example 2.1: Finite Dimensional Matrices

In this example we consider the system of ordinary differential equations where the matrix A has complex eigenvalues. To avoid these eigenvalues, we split A into two submatrices A1 and A2 that each have only real eigenvalues and compute an approximate solution based on the operator splitting: ## Show the matrix and compute eigenvalues

A=[1 2; -2 2]
[R, lambda]=eig(A);
disp('Eigenvalues of A:'), disp(diag(lambda)')

A =

1     2
-2     2

Eigenvalues of A:
1.5000 - 1.9365i   1.5000 + 1.9365i



## Split A into A1 and A2

Make a splitting so that the eigenvalues of the two new matrices are real

A1 = [0.5 0; -2 1]
[R,lambda1]=eig(A1); disp('Eigenvalues of A1:'), disp(diag(lambda1)')

A2 = A-A1
[R,lambda2]=eig(A2); disp('Eigenvalues of A2:'), disp(diag(lambda2)')

A1 =

0.5000         0
-2.0000    1.0000

Eigenvalues of A1:
1.0000    0.5000

A2 =

0.5000    2.0000
0    1.0000

Eigenvalues of A2:
0.5000    1.0000



## Investigate accuracy of the splitting

We compute the approximate solution at time T=pi using N splitting steps and report the splitting error as a function of N

T=pi;
N=250;
u0=[1 1]';
uexact=expm(-T*A)*u0;
error=zeros(N,1);
for nstep=1:N,
usplit=u0;
dt=T/nstep;
Astep=expm(-dt*A2)*expm(-dt*A1);
for i=1:nstep;
usplit=Astep*usplit;
end;
error(nstep)=sum(abs(uexact-usplit));
end;

semilogy(error,'o');
title('Error as a function of number of substeps'); axis tight;
ylabel('log(error)','Rotation',0); xlabel('number of steps'); 