1

I resolve the following equation :

enter image description here

To solve it, I would like to use matrix exponential :

enter image description here

I thought about 3 ways to do it :

  1. I could have missed it but Oj-Algo could have a simple way to compute exp(A) (I did not find it in MatrixStore javadoc)
  2. I get matrix D and V from EigenValue methods ([A] = [V][D][V]-1) and then I compute enter image description here Then the question that comes first is how I apply x->exp(x*t) function to all diagonal elements of D ?
  3. Last idea is basically the same as 2. but I previously store the scalar-matrix product in a new matrix ([X] = [D]*(-t)) and then I compute :

enter image description here

Can you help me find the best way/methods/class I should use ? Thank you

NB : This question is a "follow up question" : initial question

EDIT : This is what i have tried for now, Is it the best way to do it ? : 

import static org.ojalgo.function.PrimitiveFunction.EXP;


public class SolveDifferentialEquationTest
{

    private static final PhysicalStore.Factory<Double, PrimitiveDenseStore> matrixFactory = PrimitiveDenseStore.FACTORY;

    public static void main(String[] args)
    {
        SparseStore<Double> matrix;
        final PhysicalStore<Double> diagMatrix;
        final PhysicalStore<Double> eigenVectorMatrix;
        final PhysicalStore<Double> inverseEigenVectorMatrix;
        final Eigenvalue<Double> eigenvalue;
        final int time = 100;
        PhysicalStore<Double> initialVector;
        final PhysicalStore<Double> finalVector;

        int dim = 2000;
        matrix = SparseStore.PRIMITIVE.make(dim, dim);
        initialVector = matrixFactory.makeZero(dim,1);

        // fill matrix and initialVector
        //...

        //Decompose matrix

        eigenvalue = Eigenvalue.PRIMITIVE.make(matrix);
        eigenvalue.decompose(matrix);
        diagMatrix = eigenvalue.getD().copy();
        eigenVectorMatrix = eigenvalue.getV().copy();
        InverterTask<Double> inverter = InverterTask.PRIMITIVE.make(eigenVectorMatrix);
        try {
            inverseEigenVectorMatrix = inverter.invert(eigenVectorMatrix).copy();
        } catch (RecoverableCondition e) {
            throw new RuntimeException(e);
        }

        // Construct exp(Dt)
        diagMatrix.multiply(time);
        diagMatrix.modifyDiagonal(EXP);

        // Compute
        finalVector = inverseEigenVectorMatrix.multiply(diagMatrix)
                                              .multiply(eigenVectorMatrix)
                                              .multiply(initialVector)
                                              .copy();
    }
}
Johann MARTINET
  • 94
  • 1
  • 10
  • The exponential of a diagonal matrix is rhe diagonal matrix with the exponentials of the diagonal elements on the diagonal, eg exp( [d1 0 ; 0 d2]) = [ exp(d1) 0 ; 0 exp(d2)] – dmuir Apr 17 '18 at 08:25
  • Yes exactly, this is why I am able to compute D and exp(D) in 2. and 3. – Johann MARTINET Apr 17 '18 at 08:31
  • evd.getEigenvalues(); Perhaps it's easier to extract the eigenvalues as an array/list. If you want to modify the "D" matrix you have to copy it first as MatrixStore instances are immutable. – apete Apr 18 '18 at 15:31
  • Thank you for your answer, I edited my question with a possible solution. Can you tell me if it is the good way to use ojAlgo. Thanks – Johann MARTINET Apr 27 '18 at 09:34
  • Don't copy everything all the time - why do you do that. If the original matrix is symmetric (v is orthogonal) then you can replace the inverse with a simple transpose. – apete May 03 '18 at 11:18
  • The original matrix is not symmetric. For the copy, I will remove vector and V copy, for D I think I should keep it as you suggested in previous comments, thank you! This step is performance critical in my program, I must optimize it as much as possible. Is there a way to "give" more threads to OjAlgo (these computations will run on ~ 100 cores) ? – Johann MARTINET May 03 '18 at 14:17

0 Answers0