kill(all)$

maxima : orthogonal collocation BVP :ODE ,iso TRAM, vesion 1.0
peter.vlasschaert@gmail.com,27/11/2020

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

load(lrats) $ ratprint:false $

number of equations (without BC =( B)oundary (C)ondition)
=============================================

N:2;

polynoom & derivatives from the polynoom
=================================

≥  step 1: definitions

powerdisp:true$

→ polynoom (with x[i])

pd0[i]:=y[i]=ev(sum((x[i])^(j-1)*d[j],j,1,N+2),sum);

→ polynoom (x[1])

pd0[1];

→ first derivative : polynoom (with respect to x[i])

pd1[i]:=yd[i]=rhs(diff(pd0[i],x[i]));

→ first derivative : x[1]

pd1[1];

→ second derivative : polynoom (with respect to x[i])

pd2[i]:=yd2[i]=rhs(diff(pd1[i],x[i]));

→ second derivative : x[1]

pd2[1];

≥  step 2: y  values  ( = for collocation points), matrix

ya:makelist(pd0[i],i,1,N+2);

use :pd0[1]

→ find d[i]

yb1:subst([x[1]=1],rhs(pd0[1]));

yb2:listofvars(yb1);

yb3:apply(matrix,[yb2]);

→ find y[i]

yc1:makelist(lhs(ya[i]),i,1,length(yb2));

yc2:apply(matrix,[yc1]);

→ Q : matrix

Q:(-1)*coefmatrix(ya,yb2);

→ y[i], i =1,2..N+2  ( y=Q.d)

M1:transpose(yc1)=Q.transpose(yb2);

≥  step 3: yd (= first derivative of x[i])  values  ( = for collocation points), matrix

yda1:makelist(pd1[i],i,1,N+2);

→ C : matrix

C:(-1)*coefmatrix(yda1,yb2);

yda2:makelist(lhs(yda1[i]),i,1,length(yb2));

yda3:apply(matrix,[yda2]);

→ yd[i], i =1,2..N+2  ( yd=C.d)

M2:transpose(yda3)=C.transpose(yb2);

≥  step 4: yd2 (= second derivative of x[i])  values  ( = for collocation points), matrix

ydda1:makelist(pd2[i],i,1,N+2);

→ D: matrix

D:(-1)*coefmatrix(ydda1,yb2);

ydda2:makelist(lhs(ydda1[i]),i,1,length(yb2));

ydda3:apply(matrix,[ydda2]);

→ ydd[i], i =1,2..N+2  ( ydd=D.d)

M3:transpose(ydda3)=D.transpose(yb2);

≥  step 4: d[i] => unknow ( → d=Q^(-1)*y), matrix A,B

math:(1)
====

→ find :yd   ( → yd  =C.d = C.Q^(-1) =A.y
     A[i,j], i = internal collocations points = 2..N+1

→ find :yd2 ( → ydd=D.d = D.Q^(-1) =B.y
      B[i,j], i = internal collocations points = 2..N+1
-------------------------------------------------------------------------------------------
domain : BC1 (i=1) < <internal collocations points < <BC2 (i=N+2)
-------------------------------------------------------------------------------------------

values  : x[i] → "roots : shifted legendre polynomials ,see :website explain this roots"

A:C.invert(Q)$    /* square matrix */
B:D.invert(Q)$    /* square matrix */

A[2,1]$ /* example find : matrix element A[i,j]*/

math:(2)
====
x[i]:collocations points ' in domain: solve ISO TRAM problem '

load packages
============
see :http://maxima.sourceforge.net/docs/manual/de/maxima_65.html
see :https://en.wikipedia.org/wiki/Legendre_polynomials
ref:
new: 03/09/2019:ortho collocation method :  shifted legendre polynomials,  'version 1.0',part 1
download  (Maxima,code)       :orthocolshiftedlegendrepolynomialszeroto1version1finalpart1
download  (Maxima,png)         :orthocolshiftedlegendrepolynomialszeroto1version1finalpart1outp
new: 17/09/2019:ortho collocation method :  shifted legendre polynomials,  'version 1.0',part 2
download  (Maxima,code)       :orthocolshiftedlegendredifferentialeqversion1finalpart2
download  (Maxima,png)         :orthocolshiftedlegendredifferentialeqversion1finalpart2outp
new: 17/09/2019:ortho collocation method :  shifted legendre polynomials (hypergeometric :2F1),  'version 1.0',part3
download  (Maxima,code)       :orthocolshiftedlegendredifferentialeqversion1finalpart3
download  (Maxima,png)         :orthocolshiftedlegendredifferentialeqversion1finalpart3outp
new: 09/11/2020:maxima : appendix : 'ortho collocation method :sum,legendre',version 1.0
download  (Maxima,code)       :maximaappendixorthocollocationsumlegversion1final
download  (Maxima,png)         :maximaappendixorthocollocationsumlegversion1finaloutp

load(orthopoly)$

first : legrendre polynoom
------------------------------------

z1a:legendre_p (3, x);

z2a:factorsum(ratsimp(z1a));

second :shifted legrendre polynoom
--------------------------------------------------

z1b:legshifted(w,z):=ratsimp(legendre_p (w,2*z-1));

z1b1:legshifted(N,z);

z1b2:allroots(z1b1); /* domain : 0≤ x[i] ≤ 1 */

find :x[i], i=2,3, ....N+1
------------------------------

z1c1:makelist(x[i]=rhs(z1b2[i-1]),i,2,N+1);

→ : BC 1 & BC2,see domain

z1c2:[x[1]=0,x[N+2]=1];

z1c3:flatten(cons(z1c1,z1c2)); /* exit endcons : command */

z1c3a:sort(z1c3);

→ {A,B} =>  { AA,BB} ,            ( fill in the values,x[i],i=1,...,N+2)
------------------------------------------------------------------------------------

AA:lratsubst(z1c3,A),float;
BB:lratsubst(z1c3,B),float;

powerdisp:false;

Orthogonal Collocation 
==================
ISOTHERMAL  TRAM ( ODE : BVP (BC1,BC2))
oc :(o)rthogonal (c)ollocation method

→  ODE 
------------
Pe : Peclet Number
Da :Damkohler Number

parameter:
---------------

Pe : 6$
Da : 2$

BVP: (B)oundary (V)alue (P)roblem
--------

→ ode (ordinary differential equation, nonlinear)

od1:1/Pe*'diff(Y,x,2)-'diff(Y,x)-Da*Y^2 =0;

→ BC (BC1,BC2), Boundary Conditions

→ BC1 @ x= 0 → x[1] =0

bc1:'diff(Y,x)=Pe*(Y-1);

→  BC2 @ x= 1 → x[N+2] =1

bc2:'diff(Y,x)=0;

How to implement  :{AA,BB}  into {od1,bc1,bc2}
=====================================
----------------------------------------------------------------------------------------------------
representation : 1e order differentation : AA[i,j] ; i=1,2,...,N+1," 'diff(Y,x,1)"
                            2e order differentation : BB[i,j] ; i=1,2,...,N+1," 'diff(Y,x,2)"
----------------------------------------------------------------------------------------------------
j = 1,2,..,N+2
----------------------------------------------------------------------------------------------------
Y → y[i] into "od1"

2e order differentation
-------------------------------

u2[i]:=part(od1,1,1,1)=sum(BB[i,j]*y[j],j,1,N+2);

1e order differentation
-------------------------------

u1[i]:=part(od1,1,2,1)=sum(AA[i,j]*y[j],j,1,N+2);

system : BVP(ODE) → nonlinear equations

BC : result 2, nonlinear equations 
==========================

→ uu[1],uu[N+2]

u1[1];

uu[1]:rhs(u1[1]) = subst([Y=y[1]],rhs(bc1));

u1[N+2];

uu[N+2]:rhs(u1[N+2]) = subst([Y=y[N+2]],rhs(bc2));

fpprec : 8;

BVP(ODE) uu[i],i = 2,3,..,N+1
======================

uod[i]:= 1/6*sum(BB[i,j]*y[j],j,1,N+2) - sum(AA[i,j]*y[j],j,1,N+2)=Da*(y[i])^2;

uuu:makelist(uod[k],k,2,N+1),float,expand;

uuu[1];

uuu[2];

find : system "nonlinear equations"
===========================

→ system equations : list

sol:flatten([uu[1],makelist(uuu[k],k,1,length(uuu)),uu[N+2]]);

rhs : system of equations

sq1:makelist(rhs(sol[i]),i,1,length(sol));

sq2:apply(matrix,[sq1]);

→ system equations : matrix

M4:coefmatrix(sol,yc1),float;

M:M4.transpose(yc1) = transpose(sq2); 

numerical solution: system, mnewton
--------------------------------------------------

How to use : mnewton.
-------------------------------

ref :maxima : Tram (= finite difference), with irreversible sec order reaction,version 1.0
peter.vlasschaert@gmail.com,24/11/2020

load("mnewton")$

→ initial values (when doesn't cvg take other values )

nu1:makelist(0,k,1,N+2);

→  system equations ,list

nu2:sol;

→  variable , y[i], i=1,2,...N+2

nu3:yc1;

nu:flatten(mnewton(nu2,nu3,nu1));
              

wxplot2d: result (numerical nonlinear system).
===================================

-----------------------------------------------------------------------
see :graph,check result against the finite difference
------------------------------------------------------------------------

maxima : Tram (= finite difference), with irreversible sec order reaction,version 1.0
peter.vlasschaert@gmail.com,24/11/2020

→ x[i], i=1,2,..,N+2 ( values :shifted legendre)

vt1:z1c3a;

vt1a:makelist(rhs(vt1[i]),i,1,length(vt1));

→ y[i], i=1,2,..,N+2

vt2:nu;

vt2a:makelist(rhs(vt2[i]),i,1,length(vt2));

→ graph (wxplot2d)

l:[discrete, vt1a,vt2a];

wxplot2d([l,l], [x,0,1],[y,0,1],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [title,"orthogonal collocation (oc) ,isothermal TRAM"],        
        [legend, "discrete","continue"],
        [xlabel, "x values"],
        [ylabel, "y values"],grid2d)$