kill(all)$

maxima:LU :decomposition A=L.U , version 1.0,part 2
peter.vlasschaert@gmail.com, 28/08/2020

see website:maxima:LU :decomposition A=L.U , version 1.0,part 1

load packages,option
=================

load(lrats)$ ratprint:false$

example : A1*x =b1 
===============

simple : example, "linear system of equations"

ne:3$ /*number of equations */

pa[1]:xa+xb+2*xc = .5;
pa[2]:2*xa-xb+2*xc = 1;
pa[3]:4*xa+xb+4*xc = -1.5;

p1:makelist(pa[i],i,1,ne);

p1a:listofvars(p1);

p1aa:[makelist(p1a[i],i,1,ne)];

A1:coefmatrix(p1,p1a);

A11:augcoefmatrix(p1,p1a);

b1:-col(A11,4);

p1b:apply('matrix,p1aa);

p1c:transpose(p1b);

p1c:A1.p1b;

p1d:p1c=b1;

solve pa[i], i =1,2,3
--------------------------

p1a1:algsys(p1,p1a);

p1a2:flatten(solve(p1,p1a)),float;

use : result part 1:
==============
U,L,solve = L.U

U:matrix(
        [u[1,1],    u[1,2], u[1,3]],
        [0, u[2,2], u[2,3]],
        [0, 0,  u[3,3]]
    );

L:matrix(
		[1,	0,	0],
		[l[2,1],	1,	0],
		[l[3,1],	l[3,2],	1]
	);

M:[l[2,1]=H[2,1]/H[1,1],l[3,1]=H[3,1]/H[1,1],
l[3,2]=(H[1,1]*H[3,2]-H[1,2]*H[3,1])/(H[1,1]*H[2,2]-H[1,2]*H[2,1]),
u[1,1]=H[1,1],u[1,2]=H[1,2],u[1,3]=H[1,3],u[2,2]=(H[1,1]*H[2,2]-H[1,2]*H[2,1])/H[1,1],
u[2,3]=(H[1,1]*H[2,3]-H[1,3]*H[2,1])/H[1,1],
u[3,3]=((H[1,1]*H[2,2]-H[1,2]*H[2,1])*H[3,3]+(H[1,3]*H[2,1]-H[1,1]*H[2,3])*H[3,2]+(H[1,2]*H[2,3]-H[1,3]*H[2,2])
*H[3,1])/(H[1,1]*H[2,2]-H[1,2]*H[2,1])];

put in values : 'from example'
======================
M => M1
L  => L1
U  => U1

for i:1 thru ne do(
    for j:1 thru ne do(
        H[i,j]:A1[i,j]
    ));

M1:ev(M,nouns);

L1:lratsubst(M1,L);
U1:lratsubst(M1,U);

A1 = L1.U1
=========
A1 => A2

A2:L1.U1;

check : LU decomposition
====================
use : A1,A2

is(A1=A2);

=>   Maxima : find LU decomposition
============================

z1:lu_factor(A1,generalring)$
[D2,L2,U2]:get_lu_factors(%);

diagonal 
------------

D2;

  =============================================================================================

General theory :solve example ' xa,xb,xc'
===============================
Ax =b =(L.U).x =L.(U.x) = b

rhs : vector 'b1 ' example 'system of equations'
----------------------------------------------------------------

p2a:b1;
p2b:b1[1];

b3:[0.5,1,-1.5];

to solve : 'system of equations:simple example

need : ones (diagonal)

p2c:makelist(l[i,i]=1,i,1,ne);

combination ones  & diagonal

→use : subst,lratsubst

M2:append(p2c,M1);

→use : direct evaluation

M2a:sort([ll[1,1]:1,ll[2,2]:1,ll[3,3]:1,ll[2,1]:2,ll[3,1]:4,ll[3,2]:1,
    uu[1,1]:1,uu[1,2]:1,uu[1,3]:2,uu[2,2]:-3,uu[2,3]:-2,uu[3,3]:-2]);

1e) part solve : L.y = b , y find by forward subst.
 ----------------------------------------------------------------
y[1] =b1[1]/ll[1,1]                                                   (1)
 ----------------------------------------------------------------                                       
for i = 2. to ne
  ( y[i] = 1/ll[i,i]*(b1[i] - sum(ll[i,j]*y[j],j,1,i-1)) )      (2)
---------------------------------------------------------------- 

y[1]:lratsubst(M2,b1[1]/l[1,1]);                     /*  (1) */

for i:2 thru ne do
 y[i]:1/ll[i,i]*(b3[i] - sum(ll[i,j]*y[j],j,1,i-1)) /*  (2)  */
    ;

flatten(makelist(y[i],i,1,ne));

2e) part solve : U.x = y ,x  find by backward subst..
-------------------------------------------------------------------
x[ne] =y[ne]   /uu[ne,ne]                                        (3)
 ---------------------------------------------------------------- 
for i = ne-1. to 1 step -1
  ( x[i] = 1/uu[i,i]*(y[i] - sum(uu[i,j]*x[j],j,i+1,ne) )  (4)
----------------------------------------------------------------- 

x[ne]:y[ne]/uu[ne,ne]                               /*  (3)  */;

for i:ne-1 step -1 thru 1 do
 x[i]:1/uu[i,i]*(y[i] - sum(uu[i,j]*x[j],j,i+1,ne)) /*  (4)  */    ;

3e) result :example: x[i],i = 1,..,ne
-----------------------------------------------

res:flatten(makelist(x[k],k,1,ne));

 ====================================================================================

=>  Maxima : solve system :example
============================
use : z1 & rhs system of equations

z2:lu_backsub(z1,transpose(b3));

z3:transpose(z2);

z4:z3[1];

z5:z4[1];

check :result
==========
? p1a2 & res & z4

result : solve ,algsys  &  res  &  z4 => true 

c1:makelist(rhs(p1a2[i]),i,1,ne);