With Matlab, I have a matrix solution to find from 2 matricial equations (size matrix is 7x7). Here the 2 equations to solve with "a" and "b" are the unknow matrices and where F1, F2, P1 and P2, D, D2, D are known. Solving "a" and "b" would allow me to build a new matrix P = a . P1 + b . P2 . (remark : D matrix is equal to : D = a.a.D1 + b.b.D2 with D1and D2 diagonal matrices) :
a.a + a.P1.b.P2^T + b.P2.a.P1^T + b.b - Id = 0(equation 1)F1.a.P1 + F1.b.P2 + F2.a.P1 + F2.b.P2 − (a.P1 + b.P2).D = 0(equation 2)
Question 1) Is there a solver in Matlab that would allow to directly find matrices "a" and "b" ?
Question 2) If there is no directly solving, I think that I have to solve each solution a(i,j) and b(i,j) by using a double loop on i and j indices and call fsolve in inner loop.
The problem is that I have to respect the order in matricial product and it is difficult to not make errors when I implement the code.
Here the 2 equations that I must solve on "a" and "b" unknown :
equation(1)
equation(2)
So I tried to implement the solving of each matrix element, i.e a(i,j) and b(i,j)
For this, I am using a symbolic array x[7 7 2] and did :
% eigenv_sp = P1
% eigenv_xc = P2
% FISH_sp = F1
% FISH_xc = F2
% FISH_eigen_sp = D1
% FISH_eigen_xc = D2
% Fisher matrices
F1 = FISH_sp;
F2 = FISH_xc;
% Transposed matrix
transp_eigenv_sp = eigenv_sp'
transp_eigenv_xc = eigenv_xc'
% Symbolic variables
x = sym([7 7 2]);
aP1 = sym([7 7]);
aP1T = sym([7 7]);
bP2 = sym([7 7]);
bP2T = sym([7 7]);
d1 = sym([7 7]);
d2 = sym([7 7]);
d = sym([7 7]);
ad1 = sym([7 7]);
bd2 = sym([7 7]);
for k=1:7
for l=1:7
aP1(k,l) = sum(x(k,1:7,1).*eigenv_sp(1:7,l));
bP2(k,l) = sum(x(k,1:7,2).*eigenv_xc(1:7,l));
aP1T(k,l) = sum(x(k,1:7,1).*transp_egeinv_sp(1:7,l));
bP2T(k,l) = sum(x(k,1:7,2).*transp_egeinv_xc(1:7,l));
a2(k,l) = sum(x(k,1:7,1).*x(1:7,l,1));
b2(k,l) = sum(x(k,1:7,2).*x(1:7,l,2));
d1(k,l) = FkSH_ekgen_sp(k,l);
d2(k,l) = FkSH_ekgen_xc(k,l);
d(k,l) = d1(k,l) + d2(k,l);
ad1(k,l) = sum(x(k,1:7,1).d1(1:7,l));
bd2(k,l) = sum(x(k,1:7,2).d2(1:7,l));
end
end
% Function that represents `(a,b)(i,j)` unknown element represented by `(y(i,j,1),y(i,j,2))`
myfun=@(x,i,j) [
% First equation
sum(x(i,1:7,1).*x(1:7,j)) + sum(aP1(i,1:7).*bP2T(1:7,j)) + sum(bP2(i,1:7).*aP1T(1:7,j)) + sum(x(i,1:7).*x(1:7,2)) - eq(i,j);
% second equation
sum(F1(i,1:7).*aP1(1:7,j)) + sum(F1(i,1:7).*bP2(1:7,j)) + sum(F2(i,1:7).*aP1(1:7,j)) + sum(F2(i,1:7).*bP2(1:7,j)) - ...
sum(aP1(i,1:7).*ad1(1:7,j)) + sum(bP2(i,1:7).*ad1(1:7,j)) + sum(adP1(i,1:7).*dP2(1:7,j)) + sum(bP2(i,1:7).*dP2(1:7,j));
]
% Solution of system of non linear equations with loop
y = zeros(7, 7, 2);
for i=1:7
for j=1:7
a0 = 1e7;
b0 = 1e7;
x0 = [ a0 b0];
y(i,j,:) = fsolve(@(x)myfun(x,i,j),x0)
end
end
But at the execution, I have the following error :
Index exceeds matrix dimensions.
Error in sym/subsref (line 814)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in compute_solving_Matricial_Equations_dev (line 60)
aP1(k,l) = sum(x(k,1:7,1).*eigenv_sp(1:7,l));
I don't understand where is the issue, I only iterate on index 1:7, if someone could see what's wrong ?
Update 1
I wonder if I have the right to use for a given index i the symbolic expression x(1:7,i,1) whereas this array is not known numerically but defined only from a symbolic point of view.
Update 2
Maybe I have a clue.
If I do an expand of the symbolic array x (of size (7,7,2)), I get :
>> expand(sum(x(1,1:7).*x(1:7,1)))
Index exceeds matrix dimensions.
Error in sym/subsref (line 814)
R_tilde = builtin('subsref',L_tilde,Idx);
So, the only thing to done I think is to explicitly write the 7x7 equations
with 7x7 = 49 couples unknown x(i,j,1) and x(i,j,2) which corresponds to a(i,j) and b(i,j) matrices.
Update 3
I think the key point is to put into function myfun the 2 matrical equations but with taking into account of the 49 (7x7) couples (a(i,j) and b(i,j)) unknown that I defined previously as symbolic variables with array x[7 7 2].
And second point, I must to "make understand" to fsolve that we are searching for each iteration (double loop on (i,j)) the value of a(i,j) and b(i,j) given the whole 49 equations. Maybe it would be necessary to expand for example the expression of x(i,1:7) or x(1:7,j) in myfun function but I don't know how to proceed.
