Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.

Matrix exponential

Please login with a confirmed email address before reporting spam

Dear folks,

I have a system of ODE of the form y_dot = A(y) * y,
where y_dot is time derivative of vector y and A(y) a tensor (n by n Matrix form).
For numerical approximation, I need to calculate exponential of matrix A, which is exp(del_t * A)
where del_t is time step.
Do you have any idea how to implement this using the PDE module? or is there any other method to do this?

Thank you.
Wonseok Yoon

4 Replies Last Post 2013年6月30日 GMT-4 14:45
Jim Freels mechanical side of nuclear engineering, multiphysics analysis, COMSOL specialist

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago 2010年5月12日 GMT-4 21:57
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.
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: http://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.

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago 2010年5月16日 GMT-4 21:01
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.



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.

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago 2010年5月16日 GMT-4 21:09
Thank both of you for replies. I would try Matlab first.
Best,
Thank both of you for replies. I would try Matlab first. Best,

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago 2013年6月30日 GMT-4 14:45
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
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

Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.