Discussion Forum

In 1977-1979, I was a full-time graduate student working on a Master's in Nuclear Engineering at UT-Knoxville. We regularly used a code called "MATEXP" which solved linear ODEs using the matrix exponential method you describe. The code was written by a colleague still working today at ORNL. I don't know how you might be able to get the code today, but I did a quick search at RSICC and found it is used as a subroutine in another code: www-rsicc.ornl.gov/codes/ccc/ccc4/ccc-439.html

I am sure you could translate the essence of this code or start from scratch and write a MATEXP kernel in MATLAB and then interface that with COMSOL. You might also be able to write a MATEXP equivalent with COMSOL primitives to arrive at the equivalent.

But, this is all unnecessary because you can also solve sets of linear ODEs using other solver tools built into COMSOL.

If it's not possible to do this from the Comsol side you could probably pass the info to
Matlab, do the computation from there, and then return the results back to Comsol.

The mathematical procedure for computing any non-trivial matrix function is first to
diagonalize the matrix into its eigenvectors and eigenvalues, and then apply the identity
from Linear Algegra that says that for a nonsingular square matrix M, any function
of M may be expanded as

f(M) = X * f(S) * X'

where X and S are related to M by the orthogonal similarity transform

M = X * S * X'

Here, X is the square matrix of the column eigenvectors of M,
X' is the adjoint (complex conjugate transpose) of X, and S
is the diagonal matrix of the corresponding eigenvalues.

In your expression

exp(del_t * A)

what you would do is first diagonalize A according to

A = X * S * X'

using the Matlab eig(...) function to find X and S, then compute

exp(del_t * A) = X * exp(del_t * S) * X'

Since S is diagonal, the exp(...) term on the RHS is also diagonal
and is easily computed as just the exponential of the scalar values
along the diagonal. Completing the calculation then involves only
two additional matrix multiplications with X and X'.

If A is very large and you only need a few of the largest eigenvalues, you
can use Singular Value Decomposition (Matlab function svd(...)) to only get
the most significant elements of the matrix, as determined from the
magnitudes of the eigenvalues. The decomposition is similar to
orthogonal similarity diagonalization, but takes the form

M = U * S * V'

where U and V are not square, and where the matrix S is diagonal.
The corresponding expression for f(M) involving U and V is a little more
complicated than given above, so what you do is use svd(...) to identify
the largest eigenvalues, then use the Matlab function eigs(...) to return
only the elements corresponding to these eigenvalues. Then you can
apply the identity above directly.

You'll probably want to consult Comsol documentation about how to pass
info back and forth from Matlab, and how to invoke Matlab functions, as well as
the Matlab documentation about eig(...) and svd(...). There are things like
stability and condition numbers, etc., that you will likely have to take into
consideration.

Thank both of you for replies. I would try Matlab first.
Best,

Dear all,

I would like to define an exponential matrix in function of comsol state field variable... for example..

%Definition of matrix function
A=[Tx,0,0;0,Ty,0;0,0,Tz];

% Computation of matrix exponential
B=expm(A);

where,
T as temperature
Tx as \partial T/ \partial x
Ty as \partial T/ \partial y
Tz as \partial T/ \partial z

Matlab can perfectly compute the matrix exponential using expm()... but it uses the matrix form... and the matrix entries are the comsol expressions for me ...

Any idea or suggestion??

Your quick and kind reply is hightly appreciated
Best regards
Hamidréza