The section introduces the functions used in linear algebra and solving equations.

Linear Algebra #

inv #

Description
Computes the inverse of a matrix.
Used in FDTD and FDE.

inv( A ) computes the inverse of A. The inv function, because it uses the solve function, uses a symmetric/hermitian solver if the coefficient matrix is symmetric/hermitian.

Syntax

See also
pinv

pinv #

Description
Computes the Moore-Penrose pseudoinverse.
Used in FDTD and FDE.

Syntax

See also
inv

trace #

Description
Computes the trace of a matrix.
Used in FDTD and FDE.

Syntax

See also
diag

rcond #

Description
Condition number.
Used in FDTD and FDE.

rcond computes an estimate of the condition number of the input matrix, A. rcond() uses the LAPACK routines DGECON, or ZGECON. Probably the most published way to compute the condition of a matrix is: Kcond = ||A|| * ||inv(A)||. Another method is to use the 1st and last singular values of A, Kcond=sigma(1)/sigma(n); rcond computes an ESTIMATE of the condition number without computing all of the columns of inv(A). For more information see the LAPACK User's Guide.

Syntax

det #

Description
Matrix determinant.
Used in FDTD and FDE.

det computes the determinant of the matrix argument. det uses the LAPACK functions to factor the input, and the LINPACK algorithm to calculate the determinant.

Syntax

Example

a=[1 ,   2 ,  4;
   7,   4 ,  8;
   2 ,   3 ,  9];
d=det(a) 

Result:

d =
           -30 

rank #

Description
Rank of matrix.
Used in FDTD and FDE.

Syntax

Example

a=[1 ,   2 ,  4;
   7,   4 ,  8;
   2 ,   3 ,  9];
k=rank(a) 

Result:

k =
             3  

tril #

Description
Returns the lower triangular part of A.
Used in FDTD and FDE.

Syntax

See also
triu

triu #

Description
Returns the upper triangular part of A.
Used in FDTD and FDE.

Syntax

See also
tril

symm #

Description
symm(A) is the symmetric (Hermitian) part of A, (A+A')/2. It is the nearest symmetric (Hermitian) matrix to A in both the 2- and the Frobenius norms.
Used in FDTD and FDE.

Syntax

Example

A = [1 + 2i, 2 + 5i;...
     3 + 7i, 4 + 1i ]

A_symm = symm(A)

# A_symm is a Hermitian matrix
D = A_symm - conj(A_symm.') # or D = A_symm - A_symm'

Result:

A =
              1 +           2i              2 +           5i
              3 +           7i              4 +           1i
A_symm =
              1 +           0i            2.5 -           1i
            2.5 +           1i              4 +           0i
D =
              0 +           0i              0 +           0i
              0 +           0i              0 +           0i

permute #

Description
Rearranges dimensions of N-D array.
Used in FDTD and FDE.

This function can rearrange dimensions of a matrix according to the specified order, which can be seen as the generalization of transposition. All the elements of the input Order must be unique, real, positive, integer values.
If the input Order is [k, m, n, ... ], the size of the output matrix B will be k * m * n * ... .
And the ith dimension of the output matrix B is the Order(i)th dimension of the input matrix A .

Syntax

Example

# create matrix A
A = [1, 2, 3;...
    -2, 5, -3 ];
 
# use permute to transpose matrix A
# the 1st dimension of B is the 2nd (Order(1)) dimension of A
# the 2nd dimension of B is the 1st (Order(2)) dimension of B 
Order = [2, 1]; 
B =  permute(A, Order)    

Result:

A =
             1             2             3
            -2             5            -3
B =
             1            -2
             2             5
             3            -3

sparse #

Description
Converts full (dense) storage to sparse storage.
Used in FDTD and FDE.

sparse converts its argument from a dense storage format to the sparse storage format. If the argument is already sparse, then it is condensed (any non-zeros are removed). The sparse storage format is commonly referred to as sparse row-wise storage. Only the non-zero elements of the matrix are stored in a row-wise fashion. Row-wise storage is used for several reasons:

• The matrix vector product A * x is a very common operation, efficiently performed with row-wise storage.
• Row-wise (and column-wise) storage is a very general storage scheme that works well for general non-symmetric matrices. There is a penalty to pay for storing symmetric matrices in this fashion, but it is small.

Certain operations on sparse matrices will return dense matrices. For instance, the cosine operation on a sparse matrix will create a dense matrix with ones where there used to be zeros.

Sparse matrices are printed differently than full, or dense matrices. Only the non-zero elements are printed, along with their row and column values.

Syntax

Example

a = [0, 1, 0;
     2, 0, 0;
     0, 0, 3];
s = sparse(a)

Result:

s =
 (1, 2)    	           1
 (2, 1)    	           2
 (3, 3)    	           3

See also
spconvert, full

spconvert #

Description
Converts a full column matrix to sparse storage.
Used in FDTD and FDE.

spconvert converts its argument to, or from, the sparse storage format. If the argument is a 3 (or 4) column full matrix, the argument is converted to sparse storage format. The 1st two columns are taken as the row and column indices for the elements in the third column. The rows of the input matrix do not have to be in any particular order. If there are duplicate elements (same row and column number), then they are summed.

If the argument is a sparse matrix, then it is converted to a full matrix with 3 columns. The first two columns are the row and column indices of each non-zero element, and the third column contains the corresponding element values (columns 3 and 4 if the matrix is complex).

Syntax

Example
Creates a sparse matrix of zeros with 1000 rows, and 1000 columns.

s = spconvert([ 1000, 1000, 0 ])
show(s);

Result:

s =
	[all zeros]
	nr                  :	1000
	nc                  :	1000
	n                   :	1e+06
	nnz                 :	0
	class               :	num
	type                :	real
	storage             :	sparse

See also
sparse, full

full #

Description
Converts sparse storage to full (dense) storage.
Used in FDTD and FDE.

full converts its argument from the sparse storage format to the full, or dense, storage format.

Syntax

Example

d = [1, 1, 10; 2, 4, 20; 3, 1, 12; 5, 2, 13; 1, 4, 3]
s = spconvert(d)
f = full(s)

Result:

d =
             1             1            10
             2             4            20
             3             1            12
             5             2            13
             1             4             3
s =
 (1, 1)    	          10
 (1, 4)    	           3
 (2, 4)    	          20
 (3, 1)    	          12
 (5, 2)    	          13
f =
            10             0             0             3
             0             0             0            20
            12             0             0             0
             0             0             0             0
             0            13             0             0

See also
sparse, spconvert

diag #

Description
Constructs diagonal matrix using input argument.
Used in FDTD and FDE.

If the 1st argument A is a 1xN matrix, constructs a diagonal matrix from the input. Optionally if K (scalar) is specified then create a matrix with the vector as the Kth diagonal.

K = 0, the elements of input are arranged on the main diagonal.
K < 0, the elements of input are arranged below the main diagonal.
K > 0, the elements of input are arranged above the main diagonal.

If the 1st argument is a MxN matrix, constructs a 1xN matrix from the diagonal elements of the input matrix. Optionally if K is specified return the vector from the Kth diagonal of the input matrix.

K = 0, the elements of the newly constructed vector are from the main diagonal of the input.
K < 0, the elements of the newly constructed vector are from the diagonal below the main diagonal of the input.
K > 0, the elements of the newly constructed vector are from the diagonal above the main diagonal of the input.

Syntax

Example

# create matrix A
A = [1, 2, 3 ];

# use A to construct a new matrix 
# and the elements of A are arranged on the main diagonal 
B_onMainDiag = diag(A, 0)

# use A to construct a new matrix and
# the elements of A are arranged on the first diagonal above the main diagonal
B_aboveMainDialog = diag(A, 1)

# use A to construct a new matrix and
# the elements of A are arranged on the second diagonal below the main diagonal
B_belowMainDialog = diag(A, -2)

Result:

B_onMainDiag =
             1             0             0
             0             2             0
             0             0             3
B_aboveMainDialog =
             0             1             0             0
             0             0             2             0
             0             0             0             3
             0             0             0             0
B_belowMainDialog =
             0             0             0             0             0
             0             0             0             0             0
             1             0             0             0             0
             0             2             0             0             0
             0             0             3             0             0

See also
tril, triu

factor #

Description
Factors a square matrix.
Used in FDTD and FDE.

The factor function computes the LU factorization of the input matrix A. factor returns 3 elements.

Syntax

The user can override factor's choice of solution type with the optional argument 'TYPE'.

factor utilizes the LAPACK subroutines DSYTRF, DSYCON or ZHETRF, ZHECON.
factor returns the results in the above format, so that they may be conveniently used with backsub for repetitive solutions. The user-function lu will separate the results from factor into separate L and U matrices.

Example

A = [10, -7, 0;
     -3,  2, 6;
      5, -1, 5];

[LU, Pvt, RCond] = factor(A);
[L, U, P] = lu(A);
# LU is combined from matrix L and matrix U.
L
U
LU

Result:

L =
             1             0             0
           0.5             1             0
          -0.3         -0.04             1
U =
            10            -7             0
             0           2.5             5
             0             0           6.2
LU =
            10            -7             0
           0.5           2.5             5
          -0.3         -0.04           6.2

svd #

Description
Singular Value Decomposition.
Used in FDTD and FDE.

The output is a struct containing the three afore-mentioned objects (u, sigma, v). Various forms of the right and left singular vectors can be computed, depending upon the value of the second, optional, string argument.

Syntax

Example

A = [0.96, 1.72; 2.28, 0.96];
[u,sigma,v] = svd(A)

Result:

u =
          -0.6          -0.8
          -0.8           0.6

sigma =
             3             0
             0             1

v =
          -0.8           0.6
          -0.6          -0.8

hess #

Description
Finds the Hessenberg form of a matrix.
Used in FDTD and FDE.

Produces a unitary matrix P and a Hessenberg matrix H so that:

$A = P * H * P^{'} \\ P^{'} * P = eye(size(P))$

Syntax

Example

a = [2,3,1;45,3,5;8,9,6]
[p,h] = hess(a)

p*h*(p.')

p*p.'

Result:

a =
             2             3             1
            45             3             5
             8             9             6
p =
             1             0             0
             0  -0.984562508  -0.175033335
             0  -0.175033335   0.984562508

h =
             2   -3.12872086   0.459462504
   -45.7055795    5.50454763   -5.08808042
             0   -9.08808042    3.49545237
val =
             2             3             1
            45             3             5
             8             9             6
val =
             1             0             0
             0             1             0
             0             0             1

balance #

Description
Balances a matrix for equal row and column norms.
Used in FDTD and FDE.

balance uses the LAPACK subroutines DGEBAL and ZGEBAL to balance the input matrix so that the row and column norms are approximately equal. balance returns t and ab.

Syntax

Example

a =[ 0,  0, 1, 0;0,0,0,1;11,10,0,0;10,11,0,0]
   

[t,ab] = balance(a)

inv(t)*a*t - ab

Result:

a =
             0             0             1             0
             0             0             0             1
            11            10             0             0
            10            11             0             0
t =
             1             0             0             0
             0             1             0             0
             0             0             1             0
             0             0             0             1

ab =
             0             0             1             0
             0             0             0             1
            11            10             0             0
            10            11             0             0
            
val =
             0             0             0             0
             0             0             0             0
             0             0             0             0
             0             0             0             0            

Only square matrices are allowed.

qr #

Description
QR decomposition.
Used in FDTD and FDE.

Syntax

Example

a=[1, 2, 1, 1; 1, 1, 1, 1; 11, 10, 1, 1; 10, 11, 1, 1]
[q,r] = qr(a);

q*r

Result:

a =
             1             2             1             1
             1             1             1             1
            11            10             1             1
            10            11             1             1
val =
             1             2             1             1
             1             1             1             1
            11            10             1             1
            10            11             1             1

schur #

Description
Schur decomposition.
Used in FDTD and FDE.

The schur function returns t and z, such that:

$A = z * t * z'$

Syntax

Example

a=[1, 2, 1, 1; 1, 1, 1, 1; 11, 10, 1, 1; 10, 11, 1, 1]
       
[z,t] = schur(a);

z * t * z.'

Result:

a =
             1             2             1             1
             1             1             1             1
            11            10             1             1
            10            11             1             1
val =
             1             2             1             1
             1             1             1             1
            11            10             1             1
            10            11             1             1 

chol #

Description
Cholesky factorization, uses the LAPACK subroutine DPOTRF and ZPOTRF.
Used in FDTD and FDE.

Syntax

chol() produces an upper triangular matrix $U$, where $U'* U = A$.

Example

a = [    54,           36,           78  
         36,           29,           67  
         78,           67,          161  ]

u=chol(a);

u.' * u

Result:

a =
            54            36            78
            36            29            67
            78            67           161
val =
            54            36            78
            36            29            67
            78            67           161
 

lu #

Description
Computes the L and U factors of the matrix A.
Used in FDTD and FDE.

The function lu is a user-function that utilizes the built-in function factor. factor does the factorization using the LAPACK functions DGETRF/ZGETRF, and DGECON/ZGECON.

Syntax

Example

A = [17    24     1     8    15
    23     5     7    14    16
     4     6    13    20    22
    10    12    19    21     3
    11    18    25     2     9]

[L, U, P] = lu(A);

P.' * L * U

Result:

A =
            17            24             1             8            15
            23             5             7            14            16
             4             6            13            20            22
            10            12            19            21             3
            11            18            25             2             9
 

val =
            17            24             1             8            15
            23             5             7            14            16
             4             6            13            20            22
            10            12            19            21             3
            11            18            25             2             9

eig #

Description
Computes the eigenvectors and eigenvalues of the input matrix.
Used in FDTD and FDE.

Syntax

Example

The generalized eigenvalue problem arises quite regularly in structures. From the second order differential equations describing a lumped mass system arise $M$ and $K$, coefficient matrices representing the mass and stiffness of the various physical degrees of freedom. The equations are formulated as follows:

$M * dx^2/dt^2 + K * x = F$

Which leads to the eigenvalue problem:

$K*v = w^2 * M * v$

For a two degree of freedom system we might have:

m = eye(2,2);
k = [5,1;1,5];
[vec , val] = eig(k, m);
          
# Test the solution
          
k * vec[:,1]    
val[1] * m * vec[:,1]

Result:

val =
   -2.82842712
    2.82842712
val =
   -2.82842712
    2.82842712

Show the properties of the solution.

          
# Properties of the solution

vec' * m * vec
 
          
vec' * k * vec

Result:

val =
             1             0
             0             1
val =
             4             0
             0             6

The eigenvalues and eigenvectors are sometimes obtained by converting the generalized problem into a standard eigenvalue problem (this is only feasible under certain conditions).

mnorm #

Description
Computes the matrix norm.
Used in FDTD and FDE.

Syntax

The input parameter 'TYPE' that can be used:

• M or m
  Returns max(abs( A )).

• 1, O or o
  Returns the 1-norm (default), the largest column max(sum(abs(A))).

• 2
  Returns the matrix 2-norm (largest singular value).

• I or i
  Returns the infinity-norm, the largest row max(sum(abs(A'))).

• F, f, E or e
  Returns the Frobenius norm.

LAPACK subroutines DLANGE and ZLANGE are used to compute all norms, except the 2-norm.

Note: If 'TYPE' is Inf (the output from inf(), for example), then norm will compute the Infinity norm of the matrix A.

See also
vpnorm

vpnorm #

Description
Computes the vector P norm.
Used in FDTD and FDE.

vpnorm computes the vector P-norm of V. The second argument is required, and specifies the value of P.

Syntax

Example

# Cauchy Inequality
a = rand(3, 1);
b = rand(3, 1);

D_InEql = dot(a, b) - vpnorm(a, 2)*vpnorm(b, 2)

# There exists a constant t > 0 such that d_i= t * c_i, (i = 1，2, 3), so the inequality will be an equality.
c = rand(3, 1);
d = 2 * c;
D_Eql = dot(c, d) - vpnorm(c, 2)*vpnorm(d, 2)

Result:

D_InEql =
  -0.212200821
D_Eql =
-4.4408921e-16

norm #

Description
Computes the matrix or vector norm.
Used in FDTD and FDE.

This function is a cover-function for the two builtin functions mnorm(), which computes matrix-norms, and vpnorm(), which computes vector P-norms.

Syntax

If A is a matrix, valid B arguments are:

• 1
  Returns max(sum(abs(A))).

• 2
  Returns max(svd(A)).

• inf
  Returns max(sum(abs(A'))).

• 'fro'
  Returns mnorm(A,'F').

If A is a vector(either row or column), valid B arguments are:

• 1
  Returns sum(abs(A)).

• 2
  Returns sum(abs(A)\.^2 )^(1/2).

• p
  Returns sum(abs(A)\.^p )^(1/p).

• inf
  Returns max(abs(A)).

• - inf
  Returns min(abs(A)).

• 'fro'
  Returns sqrt(sum(diag(S * S'))).

LAPACK subroutines DLANGE and ZLANGE are used to compute all norms, except the 2-norm.

Note: If B is inf (the output from inf(), for example), then norm will compute the infinity norm of the matrix A.

Example
Example 1:

X = [2 0 1;-1 1 0;-3 3 0;-3 3 0];
n = norm(X)
n = norm(X,1)
n = norm(X,2)
n = norm(X,inf)
# sparse matrix
X = diag(ones(1,5));
n = norm(X,'fro')

Result 1:

n =
    6.33730681
n =
             9
n =
    6.33730681
n =
             6
n =
    2.23606798

Example 2:

X = sin([2 0 1 -1 1 0]);
n = norm(X)
n = norm(X,1)
n = norm(X,2)
n = norm(X,5)
n = norm(X,inf)
n = norm(X,-inf)
n = norm(X,'fro')

Result 2:

n =
    1.71785973
n =
    3.43371038
n =
    1.71785973
n =
    1.13544885
n =
   0.909297427
n =
             0
n =
    1.71785973

See also
mnorm

Solving Equations #

solve #

Description
Solves linear equations.
Used in FDTD and FDE.

solve solves a system of linear equations:

$A * X = B$

A is the coefficient matrix. B is the right hand side. X is the solution.
solve uses the LAPACK subroutines DGETRF, and ZGETRF if A is general. solve uses the LAPACK subroutines DSYTRF, and ZHETRF if A is symmetric.

Syntax

The third and optional argument 'TYPE' allows the user to override the symmetry check, and force the solve to use either the general or the symmetric solver:

• If TYPE is 'g' or 'G', the general solution is used.
• If TYPE is 's' or 'S', the symmetric solution is used.

Example

A = [1, 0;1, 1]

B = [2; 3]

X = solve(A, B)

Pr = A * X

Result:

A =
             1             0
             1             1
B =
             2
             3
             
X =
             2
             1
Pr =
             2
             3

sylv #

Description
Solves the Sylvester matrix equation.
Used in FDTD and FDE.

Syntax

poly #

Description
Converts roots to polynomial.
Used in FDTD and FDE.

Syntax

Example
Example 1:

# example 1
A = [1 2 3; 4 5 6; 7 8 0];
p = poly(A)

Result:

p =
             1            -6           -72           -27

Example 2:

# example 2
A = [1 8 -10; -4 2 4; -5 2 8];
p = poly(A);  
e = eig(A);

deg = length(p);

r = e(1);
y = 0;


for k = 1 : deg
    
    y = y +  p(deg + 1 - k) * (r^(k-1));
            
end    
  
y  

Result:

y =
  -1.136868377e-12 +              0i

lyap #

Description
Solves the general form of the Lyapunov (Sylvester) equation.
Used in FDTD and FDE.

lyap solves the general form of the Lyapunov (Sylvester) equation:

$A* X + X * B = -C$

To solve the special form of the Lyapunov equation:

$A*X + X*A' = -C$

We can use this method: lyap (A,A',C).
Syntax

See also
schur, sylv

backsub #

Description
Solution of $Ax = B$ by backsubstitution.
Used in FDTD and FDE.

The backsub function computes the solution to the set of linear equations described by:

$A * X = B$

The 1st argument to backsub is a struct converted from the result of factor(A). The second argument to backsub is the matrix B. The matrix B can contain multiple right hand sides.

backsub returns a matrix X which contains the solution(s) to the aforementioned equations.

backsub utilizes the LAPACK subroutines DGETRS or ZGETRS if the input struct contains LU factors or LAPACK subroutins DSYTRS or ZHETRS if input struct contains the LDL factors.

Syntax

Example

• A is a general matrix.

A = [1, 2, 3; 4, 5, 6; 7, 8, 0];
B = [1; 2; 3];

[Lu, Pvt, Rcond] = factor(A);
F = struct('lu': Lu, 'pvt': Pvt, 'rcond': Rcond);
X = backsub(F, B)

# verify the matrix X is the solution of matrix equation A * X = B
D = A * X - B

Result:

X =
  -0.333333333
   0.666666667
             0
D =
-1.11022302e-16
             0
             0

• A is a symmetric matrix

A = [1, 2, 0; 2, 5, 6; 0, 6, 9]
B = [5; 30; 39];

[Ldl, Pvt, Rcond] = factor(A);
F = struct("ldl": Ldl, "pvt": Pvt, "rcond": Rcond);
X = backsub(F, B)

# verify the matrix X is the solution of matrix equation A * X = B
D = A * X - B

Result:

X =
             1
             2
             3
D =
             0
             0
             0

See also
factor, inv, lu, solve

ode #

Description
Integrates Ordinary Differential Equations.
Used in FDTD and FDE.

ode integrates a system of N first order ordinary differential equations of the form:

$dy(i)/dt = f(t, y(1), y(2), ..., y(N))$

y(i) given at t.

The arguments to ode are:

• rhsf
  A function that evaluates dy(i)/dt at t. The function takes two arguments and returns dy/dt. An example that generates dy/dt for Van der Pol's equation is shown below.

function xp = vdpol(t, x)
    xp = zeros(2, 1);
    xp(1) = x(1) * (1 - x(2)^2) - x(2);
    xp(2) = x(1);
end;

• ystart
  The initial values of y, y(tstart).

• tstart
  The initial value of the independent variable.

• tend
  The final value of the independent variable.

• dtout
  The output interval. The vector y will be saved at tstart, increments of tstart + dtout, and tend. If dtout is not specified, then the default is to store output at 101 values of the independent variable.

• relerr
  The relative error tolerance. Default value is 1e-6.

• abserr
  The absolute error tolerance. At each step, ode requires that:

$abs(local error) <= abs(y)*relerr + abserr$

For each component of the local error and solution vectors. The default value is 1e-6.

• uout
  Optional. A user-supplied function that computes an arbitrary output during the integration. The uout must return a row-matrix at each dtout during the integration. It is entirely up to the user what to put in the matrix. The matrix is used to build up a larger matrix of the output, with one row for each dtout. The resulting matrix is returned by ode when the integration is complete.

The Fortran source code for ode is completely explained and documented in the text, 'Computer Solution of Ordinary Differential Equations: The Initial Value Problem' by L. F. Shampine and M. K. Gordon.

Syntax

See also
backsub, filter