kill(all)$

maxima :algebra ,inverse sum of two matrix(both inverse),version 1.0
peter.vlasschaert@gmail.com, 14/10/2019

load : packages
----------------------

load(lrats)$

/*declare(A,nonscalar,B,nonscalar)*/;

/*properties(A)*/;

rule : M.M^^(-1) = 1
===============
M :matrix 
I   : general unit matrix =1,diagonal are all ones and rest of the entries are zero.
entries of matrix ( = position inside the matrix,same : 'excel sheet')
-----------------------
i: rows
j : columns

p1:(I+A)^^(-1);

rem : inverse(M)  = M^^(-1)

I ,A same dimension ,indices 
=>  number of rows = number of columns

p2:(I+A);

A,B : matrix 
----------------

rule : M.M^^(-1) = M^^(-1).M = 1
========================
. = noncommutative operator
rem : A.B ≠ B.A
rem : M :I+A  

p3:'p1.'p2 =p1.p2;

p4:'p2.'p1 =p2.p1;

rule : I = A..A^^(-1) = A^^(-1).A  = 1
         I = B..B^^(-1) = B^^(-1).B  = 1

simp :false;

p5a:I = A.A^^(-1);
p5b:I = A^^(-1).A;
p5c:I = B.B^^(-1);
p5d:I =  B^^(-1).B;

rem : 'p5e' is a list , different set
--------------------------------------------
list :[,,,,,]
set:{,,,,}
--------------------------------------------

p5e:[p5a,p5b,p5c,p5d];
p5f:rhs(p5e[1]);

rem : A,B,M  :matrix

show the identity
=============
general : lhs = rhs

(A^^(-1)+B^^(-1))^^(-1) = B.(A^^(-1)+B^^(-1)).A = A.(A^^(-1)+B^^(-1)).B

1e)  ( I+A^^(-1))  ?
-------------------------------

p6: ( I+A^^(-1)); 

a)  simp : false
--------------------

p6b:subst([p5b],p6);

p6c:p6b=A^^(-1)*(A+I );

p6d:ev((p6c)^^(-1),simp:true);

p6d1:args(rhs(p6d));

p6d2:args(p6d1[1]);

rule : (A.B)^^(-1) = B^(-1).A^(-1)
========================
suppose : X,Y :matrix

p6e:(X.Y)^^(-1) = Y^^(-1).X^^(-1);
p6e1: X =p6d2[1];
p6e2 :Y =p6d2[2];

p6e3:lratsubst([p6e1,p6e2],p6e);

p6e4:lhs(p6e3)  = rhs(p6e3);

b) result 

p6e5:lhs(p6e4) = ev(rhs(p6e4),simp:true);

rem : rule not known.

rule : (A^^(-1))^^(-1)  = A
===================
use : rule on rhs of 'p6e4'

p6f:(X^^(-1))^^(-1)  = X;

p6f1:subst([X=A],p6f);

p6f2:ev(lratsubst([p6f1],rhs(p6e4)),simp:true);

general : evaluation for simplification.
---------------------------------------------------

simp;

result : identity
--------------------

p7:(p6)^^(-1) = p6f2;

rem :  I+A = A+I
**:rule : A+B = B+ A
=====================================
rem : A,B both matrix or vector, same dimension
in rows and columns

next thing to show :

p8a:(A^^(-1)+B^^(-1))^^(-1) = B.(A^^(-1)+B^^(-1)).A; 
P8b:(A^^(-1)+B^^(-1))^^(-1) = A.(A^^(-1)+B^^(-1)).B;

because :rule of sum above **
B => A , A=> B , result : p8b

lemma 1: from above.
----------------------------

p9:(I+X^^(-1))^^(-1)  =subst([A=X],rhs(p6e5));

rem : X(general) ,I =unit matrix

1e) start :  lhs 'p8a'
==============

p9a:lhs(p8a);

p9a1:args(p9a);

first       element  : term
second element : power

first : term
========

p9a2:p9a1[1] = part(p9a1[1],1) +I* part(p9a1[1],2) ;

use : p5b

p9a3:lhs(p9a2) = subst([p5b],rhs(p9a2));

rule : X+(X.P).Q  = X.(I+P.Q)
rem :X,P,Q : matrix
rem : P.Q ≠ Q.P

use :rule

p9a4:X+(X.P).Q  = X.(I+P.Q);

use :rhs(p9a3)

p9a5a:[X=A^^(-1),P=A,Q=B^^(-1)];
p9a5b:lratsubst(p9a5a,p9a4);

p9a6:lhs(p9a2) = rhs(p9a5b);

second : power
===========

p9a7:(lhs(p9a6))^^(-1) =ev( (rhs(p9a6))^^(-1),simp:true) ;

use :rule 
rule : (X.Y)^^(-1) =Y^(-1).X^(-1)

p9a8a:[X=A^^(-1),Y=(I+A . B^^(-1))];
p9a8b:lratsubst(p9a8a,p6e);

use :rule 
rule : (X^^(-1))^^(-1)  = X

p9a9a:[X=A];
p9a9b:lratsubst(p9a9a,p6f);

1e part
-----------

p9a9c:subst([p9a9b],p9a8b);

use :rule 
rule : p5c

p9a10:rhs(p9a9c) = subst([p5c],rhs(p9a9c) );

p9a10a:args(rhs(p9a10));

p9a10b:args(p9a10a[1]);

rule :P.Z+Q.Z=(P+Q).Z
rem :Z,P,Q : matrix

p9a11a:[P=B,Q=A,Z=B^^(-1)];
p9a11b:P.Z+Q.Z=(P+Q).Z;

p9a11c:lratsubst(p9a11a,p9a11b);

2e part
-----------

p9a12:part(rhs(p9a10),1) = (rhs(p9a11c))^^(-1);

use :rule 
rule : (X.Y)^^(-1) =Y^^(-1).X^^(-1)

p9a13a:[X=A+B,Y=B^^(-1)];
p9a13b: lratsubst(p9a13a,p6e);

use :rule 
rule : (X^^(-1))^^(-1)  = X

p9a14a:[X=B];
p9a14b: lratsubst(p9a14a,p6f);

p9a14c:rhs(p9a14b).part(rhs(p9a13b),2);

rem : A+B =B+A

3e part: final result
-------------------------

p9a15a:lhs(p8a)  = p9a14c.part(rhs(p9a10),2);

rem : see above: other version

p9a15b:P8b;