Search code examples
matlabmatrixequationdifferential-equations

Solving a system of matrix equations using MATLAB?


I have a system of three equations that I'd like to solve via MATLAB, and I'm a bit confused on how to go about doing it.

I have three equations:

A = R*P1
B = R*P2
C = R*P3

A, B, C, and P1, P2, and P3 are 3x1 matrices, while R is a 3x3 matrix. R is the same for all three equations.

I need to find R, and am given A, B, C, and the P's.

I wanted to use fsolve, but it seems that fsolve doesn't work when the variables are matrices. What is an alternative method that you would recommend using?

Just making up some numbers to work with:

P1 = [1;1;1];
P2 = [2;3;4];
P3 = [5;4;3];

R = [2 4 5; 1 5 4; 1 2 3];

Which would mean that:

A = [11;10;6];
B = [36;33;20];
C = [41;37;22];

Solution

  • If A, B, C, P1, P2, P3 are all numerical, why don't you simply use the ldivide or \ operator? This will allow you to solve your linear system directly. I can see that you have the following relationships:

    R*P1 = A
    R*P2 = B
    R*P3 = C
    

    You can see that each matrix equation yields three constraints. What you can do is create one system that encapsulates all matrix equations together thus yielding 9 constraints. As such, you'll need to reformulate this to be able to solve for your coefficients in your matrix R differently. To do this, we need to reshape your matrix R such that it becomes a 9-element vector. In other words, we can then reformulate your system like so:

    [P1 0 0 0 0 0 0]    [R1]   [     ]
    [0 0 0 P1 0 0 0]    [R2]   [  A  ]
    [0 0 0 0 0 0 P1]    [R3]   [     ]
    [P2 0 0 0 0 0 0]    [R4]   [     ]
    [0 0 0 P2 0 0 0]  * [R5] = [  B  ]
    [0 0 0 0 0 0 P2]    [R6]   [     ]    
    [P3 0 0 0 0 0 0]    [R7]   [     ]
    [0 0 0 P3 0 0 0]    [R8]   [  C  ]
    [0 0 0 0 0 0 P3]    [R9]   [     ]
    
         P           *   R   =    D
    

    You'll see that we have a 9 x 9 matrix called P, our matrix R which is reshaped into a vector so that we can solve for the coefficients and D are A,B,C concatenated together into one vector. R1 to R9 are the coefficients of the matrix R that are read from left to right and top to bottom.

    Therefore, to find your coefficients in your matrix, simply do:

    R = P^{-1}*D
    

    As such, simply construct the matrix P and vector D like so:

    P = [P1.' zeros(1,6); zeros(1,3) P1.' zeros(1,3); zeros(1,6) P1.'; ...
    P2.' zeros(1,6); zeros(1,3) P2.' zeros(1,3); zeros(1,6) P2.'; ...
    P3.' zeros(1,6); zeros(1,3) P3.' zeros(1,3); zeros(1,6) P3.'];
    
    D = [A; B; C];
    

    Now, simply solve for R and reshape it back into a 3 x 3 matrix. Therefore:

    R = P \ D;
    R = reshape(R, 3, 3).';
    

    reshape turns our vector into a 3 x 3 matrix, but it constructs the matrix in column-major format, so you need to transpose the result after you call reshape. With your example, this is what we get. I construct P1, P2, P3, A, B, C, then use the code I had from before:

    P1 = [1;1;1];
    P2 = [2;3;4];
    P3 = [5;4;3];
    
    A = [11;10;6];
    B = [36;33;20];
    C = [41;37;22];
    
    P = [P1.' zeros(1,6); zeros(1,3) P1.' zeros(1,3); zeros(1,6) P1.'; ...
    P2.' zeros(1,6); zeros(1,3) P2.' zeros(1,3); zeros(1,6) P2.'; ...
    P3.' zeros(1,6); zeros(1,3) P3.' zeros(1,3); zeros(1,6) P3.'];
    
    D = [A; B; C];
    
    R = P \ D;
    R = reshape(R, 3, 3).';
    

    To verify that R is correct, do:

    A1 = R*P1;
    B1 = R*P2;
    C1 = R*P3;
    

    We get for each:

    A1 =
    
        11
        10
         6
    
    B1 =
    
        36
        33
        20
    
    C1 =
    
        41
        37
        22
    

    This matches up with your example. However, note that you may get a warning that R is ill-conditioned. This is because you may not have enough constraints to properly find a unique inverse. You may have to add more constraints to get a unique inverse, but if you can't, then use this with caution.