kill(all)$

appendix :maxima:fvm:tridiagonal system (algo:Thomas),version 1.0
peter.vlasschaert@gmail.com,21/10/2022

info : solve tridiagonal system of equations
-------------------------------------------------------

fvm = (f)inite (v)olume (m)ethod
rem :algo:Thomas =TDMA 
TDMA = Tri Diagonal (Ma)trix

used : one direction flow ( example : x-direction)

load packages
===========

load(lrats)$

option
====

ratprint:false$

parameters
=========

st:1;
n:5;/*number of equations */

1e) general tridiagonal matrix structure
============================

→ begin equation :

p1:eq[st]:=A[st]*x[st-1]+B[st]*x[st]+C[ist]*x[st+1]=D[st]$

→ general equation :

for k:2 thru n-1 do(
    p2:eq[k]:=A[k]*x[k-1]+B[k]*x[k]+C[k]*x[k+1]=D[k]
)$

→ end equation :

p3:eq[n]:=A[n]*x[n-1]+B[n]*x[n]+C[n]*x[n+1]=D[n]$

a) list :equations
----------------------

leq:makelist(eq[ka],ka,1,n);

correction : x[0]=0,x[n+1]=0

leq:lratsubst([x[0]=0,x[n+1]=0],leq);

b) list : variables
---------------------

lvar:makelist(x[k1],k1,1,n);

c) list : rhs MA
------------------

tridiagonal system 
=============

MA:coefmatrix(leq,lvar);

The coefficient matrix(=MA) which has nonzero elements on its
principal diagonal ( = superdiagonal) and two (=subdiagonal ) with
zero elements elsewhere.

MB:augcoefmatrix(leq,lvar);

algo : Thomas method ( = recurrence formula)
==================================

p[1]:C[1]/B[1];
q[1]:D[1]/B[1];

for j:2 thru n-1 do(
p[j]:C[j]/(B[j]-A[j]*p[j-1])
)$

for j:2 thru n do(
q[j]:(D[j]-A[j]*q[j-1])/(B[j]-A[j]*p[j-1])
)$

p's and q's  :computed ,'backwards' solved 
===============================
use : step -1
→  (x[n],x[n-1],x[n-1],...x[2],x[1])

x[n]:q[n];

for k:n thru 1 step -1
  do(
 if k= n then 
    x[n]:q[n]  
  else 
    x[k]:q[k]-p[k]*x[k+1]
)$

x[1];

example :
=======

leq:["system of equations ,n by n :","","","",""];

eq1[1] :4*y[1]+y[2]=1.2;
eq1[2] :y[1]-4*y[2]+y[3]=-.31;
eq1[3]:3*y[2]+4*y[3]-y[4] =1.33;
eq1[4]:-y[3]-y[4]+4*y[5] =0.85;
eq1[5]:4*y[4]-y[5] =0.59;

leq1:makelist(eq1[k2],k2,1,5);

leq2:makelist([eq1[k2]],k2,1,5);

TableForm :equations
================

lqe3:transpose(matrix(leq,leq1));

rem : leq & leq1 must be of the same dimension

table_form(lqe3)$

levar1:makelist(y[k2],k2,1,5);

ME:coefmatrix(leq1,levar1);

p4a:maplist("=",MA,ME);

length(p4a);

p4a[1];

lhs(p4a[1]);

lhs(p4a[1])[1];

rhs(p4a[1]);

rhs(p4a[1])[1];

p4a[2];

ME1:augcoefmatrix(leq1,levar1);

→ rhs :eq1[k], k:1..n

p4b:makelist(D[i]:-ME1[i,n+1],i,1,n);

(-) sign : result of construction,augmented matrix
===================================
reason : eq1[j],j:1..n , so reduced  rhs = 0
---------------------------------------------------------------

D[1];

→ lhs :eq1[k], k:1..n

    
tes[ka]:=makelist( lhs(p4a[ka])[kb]=rhs(p4a[ka])[kb],kb,1,n);   

tes[1];

p4c:makelist(tes[kc],kc,1,n);

p4c1:flatten(p4c);

1e) removes doubles ,trples,.......... from a list
=================================
→  list to set  : setify

setify(%);

2e) set back to list
=============
→  set to list :listify

listify(%);

% , see result before.

resex1:lratsubst(%,x[1]),float;

solution :example,'ME1'
=================

list : y[i]=x[i], i:1..5
=============

totres:makelist(y[kc]=lratsubst(p4c1,x[kc]),kc,1,n),float;

why : because problem  with  :

totresfinal:ev(totres,nouns);

check : eq1[i],i =1..n, solution
=====================

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

lratsubst(totresfinal,eq1[1]),float;

→  'all solution check for all equations'

makelist(lratsubst(totresfinal,eq1[kd]),kd,1,n),float;

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