kill(all)$

maxima : algo TDMA ( = Thomas algo);version 1.0
peter.vlasschaert@gmail.com,21/09/2020

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

(load(lrats),ratprint:false) $

TDMA = solve tridiagonal system of equations
====================================

1e) step : build equations, tridiagonal system.
---------------------------------------------------------------

r1:makelist(A[i]*x[i-1]+B[i]*x[i]+C[i]*x[i+1]=D[i],i,1,5);

r2:[x[0]=0,x[6]=0];

r1a:lratsubst(r2,r1);

r1b:sort(listofvars(r1a));

→   list of variables

r1c:makelist(x[i],i,1,length(r1a));

→ list of known values,rhs

r1d:makelist(D[i],i,1,length(r1a));

2e) step : tridiagonal system
------------------------------------------
use :example for distillation problems.

r3:coefmatrix(r1a,r1c);

r3a:makelist(D[i],i,1,length(r1a));

3e) algo : Thomas.
----------------------------

→  a)  C[1]*x[2]+B[1]*x[1]=D[1],use r1a[1]=x[1]

r4a:solve(r1a[1],r1c[1])[1],expand;

r4b:q[1]=lratsubst([x[2]=0],rhs(r4a));

r4c:p[1]=(-1)*ratcoef(rhs(r4a),x[2]);

matrix 'r3' : A[1] = 0, B[1] = 1, C[1]  = p[1],D[1] = q[1]

r4d:lratsubst([reverse(r4b),reverse(r4c)],r4a);

r4da:reverse(r4d);

→  b) C[2]*x[3]+B[2]*x[2]+x[1]*A[2]=D[2],use r1a[2] =x[2]

r5a:solve(r1a[2],r1c[2])[1],expand;

use :r4d

r5b:lratsubst([r4d],r5a),expand;

solve :x[2]

r5c:solve(r5b,x[2])[1],expand;

find :p[2],q[2], i=2

r5c1:p[2]=lratsubst([x[3]=0],rhs(r5c));

r5c2:q[2]=(-1)*ratcoef(rhs(r5c),x[3]);

r5c3:x[2]=q[2]-p[2]*x[3];

matrix 'r3' : A[2] = 0, B[2] = 1, C[2]  = p[2],D[2] = q[2]

→ c) p[n],q[n], general ,n,solve :x[n]

r6a:q[n]=(D[n]-q[n-1]*A[n])/(B[n]-p[n-1]*A[n]);

r6b:p[n]=C[n]/(B[n]-p[n-1]*A[n]);

rem :p[1],q[1]
-------------------

r6c1:lratsubst([n=1],r6a);

use : q[0]=0,p[0]=0

r6c1a:lratsubst([q[0]=0,p[0]=0],r6c1);

r6c2:lratsubst([n=1],r6b);

r6c2a:lratsubst([q[0]=0,p[0]=0],r6c2);

 =>end  rem 

solve:x[n]

r6c:x[n]=q[n]-p[n]*x[n+1];

matrix 'r3' : A[n] = 0, B[n] = 1, C[n]  = p[n],D[n] = q[n]

n = 1,2,....N-1, N( =n= 5)

r7:lratsubst([n=5],r6c);

→ n+1 = 6  =>  x[5+1] = 0

r7a:lratsubst([x[5+1]=0],r7);

→ d) new matrix :r3

r8:makelist([A[i]=0,B[i]=1,C[i]=p[i],D[i]=q[i]],
    i,1,length(r1c));

r8a:flatten(r8);

r4:lratsubst(r8a,r3);

r4a:r4.(r1c) =transpose(lratsubst(r8a,r3a));

→ e) solve :successive  x[n ],trough backward substtitution .
use :n → n-1,into r6c

r9:lratsubst([n=n-1],r6c); /*same result:step = -1 when use r6c */

algo :Thomas:TDMA
-----------------------------
computation : p[1],q[1]
for n=2 to N
  {computation : p[n],q[n]}
x[N] = q[N]
for n =N-1 to 1 step -1
  {x[n]=q[n]-p[n]*x[n+1]}