kill(all)$

maxima : orthogonal collocation BVP :PDE c(x,t),iso TRAM, version 1.0
peter.vlasschaert@gmail.com,01/11/2020

ref: see website 
=============
-------------------------------------------------------------------------------------------------------------------
new: 27/11/2020:maxima : orthogonal collocation BVP :ODE ,iso TRAM, vesion 1.0
youtube :Orthogonal Collocations Method for Solving ODE - BVPs and PDEs
download  (Maxima,code)       :maximaorthocollocisotramversion1final
download  (Maxima,png)         :maximaorthocollocisotramversion1finaloutp
--------------------------------------------------------------------------------------------------------------------

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

load(lrats) $ ratprint:false$ load(draw) $

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

N:4;

powerdisp:true;

→  PDE 
------------

-------------------------------------------------------
Pe : Peclet Number
Da :Damkohler Number
==============================
oc =(o)(c) = (o)rtho (c)ollocation method
------------------------------------------------------
domain : 0 ≤ x ≤ 1, x=space ,t =time

Pe:6;
Da:2;

pde (partial differential equation, time(=t) ,space(=x))
========================================
i=2,..,N+1, x[i]

pde:diff(c(x,t),t)=1/Pe*diff(c(x,t),x,2)-diff(c(x,t),x)-Da*c(x,t)^2 ;

→rem:oc we convert : pde → nonlinear ode 

→ BC (BC1,BC2), Boundary Conditions
===============================
i = 1,N+2  (begi,end)  of the domain

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

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

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

bc2:diff(c(x,t),x)=0;

→rem:oc we convert:ode →2  nonlinear equations

variables : numerical method ( =oc)

cl:makelist(c[i],i,1,N+2);

cll:makelist(sconcat(c,i),i,1,N+2);

xl:makelist(x[i],i,1,N+2);

xll:makelist(sconcat(x,i),i,1,N+2);

→ polynoom

pde0[i]:=ye[i]=ev(sum((x[i])^(j-1)*c[j],j,1,N+2),sum);

pde0[1];

→ frist derivative :polynoom

pde1[i]:=yde[i]=rhs(diff(pde0[i],x[i]));

pde1[1];

→ second derivative : polynoom

pde2[i]:=yd2e[i]=rhs(diff(pde1[i],x[i]));

pde2[1];

load(orthopoly)$

collocation points: internal ( shifted legendre)
====================

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

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

→ Q : matrix

ya:makelist(pde0[i],i,1,N+2)$

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

→ C : matrix

yad:makelist(pde1[i],i,1,N+2)$
C:coefmatrix(yad,cl);

→ D: matrix

yad2:makelist(pde2[i],i,1,N+2)$
D:coefmatrix(yad2,cl);

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

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

1e) first : input  , x[i],shifted legendre roots, Q,C,D ' =matrix '
==============================================
use :z1b2

→ internal collocations points ( i=2,...,N+1)

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

→ BC  (i=1,N+2)

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

m1:sort(flatten(cons(m1a,m1b)));

fpprec : 3;

if numer#false then numer:false else numer:true;

2e)  {C,D,Q} → {CC,DD,QQ} , conversion ,with input of values
===============================================

CC:lratsubst(m1,C),float;
DD:lratsubst(m1,D),float;
QQ:lratsubst(m1,Q),float;

3e) {A,B} 
=======

A:CC.invert(QQ);    /* square matrix */
B:DD.invert(QQ);    /* square matrix */

4e)  BC ,{c[1],c[N+2]}
================
convert :ode (BC) into (two nonlinear equations)

bc1n:sum(A[1,j]*c[j],j,1,N+2)=Pe*(c[1]-1);
bc2n:sum(A[N+2,j]*c[j],j,1,N+2)=0;

5e)  BC ,{c[1],c[N+2]}  → find :c[1],c[N+2]
==================================
solve two nonlinear equations →

bc1a:expand(solve(bc1n,c[1]));
bc2b:expand(solve(bc2n,c[N+2]));

6e) construct  : N+2-2 = N , differential equations
======================================
combine : 5e

1e)  step
=======

cpde[i]:= 1/6*sum(B[i,j]*c[j],j,1,N+2) - sum(A[i,j]*c[j],j,1,N+2)=Da*(c[i])^2;

→ i =2..N+1

c1pde:makelist(cpde[i],i,2,N+1),expand;

c2pde:makelist(lhs(c1pde[i])-rhs(c1pde[i]),i,1,N);

2e)  step , symbolic , differential equations(without BC:bc1,bc2)
==================================================

c3pde:makelist('diff(c[i],t)=c2pde[i-1],i,2,N+1);

3e)  step , use '5e) BC '
==================

c4p:flatten([bc1a,bc2b]);

→ c[6] from BC

c4p1:expand(lratsubst(bc1a,bc2b[1]));

c4p2:expand(solve(c4p1,c[N+2]));

→ c[1] from BC

c4p1a:expand(lratsubst(bc2b,bc1a[1]));
c4p2a:expand(solve(c4p1a,c[1]));

c4pp:flatten(cons(c4p2,c4p2a));

4e)  step , symbolic , differential equations(with  BC:c4p)
===========================================

c4pde:expand(lratsubst(c4pp,c3pde));

→ use for rk5f  ,

c4rpde:makelist(rhs(c4pde[i]),i,1,N);

w1:table_form(c4rpde)$

5e) step representation :differential equations d(c[i])/dt, i=2 to N+1
------------------------------------------------------------------------------------------

d1pde:transpose(apply(matrix,[(c4pde)]));

w1a:table_form(d1pde)$

7e) solve w1a with rk5f 
==================

How to use :

rkf45([k1*(C-L),k2*(32-C)+k3*(L-C),[L,C],[150,150],[t,t_start,t_end)

a)  equation rkf45
------------------------

wa:c4rpde;

wa1:makelist(cl[i]=cll[i],i,2,N+1);

wa2:subst(wa1,wa);

wx:[c[2]=s2,c[3]=s3,c[4]=s4,c[5]=s5];

wb1:expand(lratsubst(wx,wa));

b) list variables  rkf45
-------------------------------

c) initial values for rkf45
--------------------------------

wc:[0.5,0.5,0.5,0.5];

d) simulation domain
------------------------------

wd:[t,0,2.0];

e) solve with rkf45
--------------------------

load(rkf45)$
load(BDF)$

we:rkf45(-wb1,[s2,s3,s4,s5],wc,wd)$

???  :We need place wb1 → (-) wb1 , otherwise the result  : goes to → ∞, doesn't cvg 
=================================================================

--------------------------------------------------------------------------------------------------------------------
check :eigenvalues : complex values (neg :real part, when :reaction term = removed).
--------------------------------------------------------------------------------------------------------------------

see:website , not part share contrib  : maxima ( used :stiff problems, copy package their)
see:https://themaximalist.org/tag/wxdraw2d/

help(BDF2);

  => solution :pde , 'see above'
========================

we$/*length ,we */

analysis : data :' data structue :we '

a)  find : initial values t=0
===================
other values (x[2,x[3],...x[N+1])), BC are already included

------------------------------------------------------------------------------
x[i]:collocation points ,shifted legendre ,i = 1,2,...,N+2,all
------------------------------------------------------------------------------

→ begin : we

solp1:we[1];

→ initial : time (t=0)

inittime : solp1[1];

→ end : we

solp2:we[length(we)];

→ final : time (t=final time)

finaltime:solp2[1];

→ simulation time:simtime

simtime:finaltime-inittime;

b) first :graph :c[2](t),c[3](t),...,c[N+1](t), solution of  "differential equations:w1a"
===========================================================
x[i],i=2,3,....,N+1  , internal collocations points

→ ref :https://themaximalist.org/tag/wxdraw2d/
i:fuction :nth( or mth)  = {first,second ,third,fourth,fifth,sixth,seventh},see' below'.

sol_points(numsol,nth,mth):=points(map(nth,numsol),map(mth,numsol));

grp[i]:=sol_points(we,first,i);

wxdraw2d(point_type=-1,
points_joined = true,color=red,grp[second],color=blue,grp[third],
color = green,grp[fourth],color = orange ,grp[fifth],
title="pde→differential solution(=w1a),oc:method");

color:red:c[2](t),blue:c[3](t),green:c[4](t),orange:c[5](t)
rem: reason negative values ( = see above ,wb1 → -wb1),reason:graph is upsitedown.

c) second :graph :c[2](t),c[3](t),...,c[N+1](t), solution of  "differential equations:& BC"
===============================================================
i=1,2,3,....,N+1,N+2 ,all collocations points, use BC :c4pp,'see above'

→  step 1:BC

grapha:sort(c4pp);

→ step 2: c[2]

→ find :c[2], second element of the list 'we'
[0.006441711579005277,0.4608800236013914,0.5083470512014654,0.5006561897760702,0.5051818561528478]

we[2];

we[2][2];

→ step 3: find c[1] ,,values (=bc1)

kp1:cnx[i]:=lratsubst([c[2]=we[i+1][2],c[3]=we[i+1][3],c[4]=we[i+1][4], c[5]=we[i+1][5]],
    grapha[1])$            

kp2:makelist(cnx[i],i,0,length(we)-1);

length(we);

length(kp2);

→ step 4: find c[N+2],values (=bc2)

kq1:cny[i]:=lratsubst([c[2]=we[i+1][2],c[3]=we[i+1][3],c[4]=we[i+1][4], c[5]=we[i+1][5]],
    grapha[2])$
kq2:makelist(cny[i],i,0,length(we)-1);
kq3:length(kq2);

→ step 5: first element with list :we & [c[1],c[N+2]]

za1:[we[1][1],rhs(kp2[1]),we[1][2],we[1][3],we[1][3],we[1][4],we[1][5],rhs(kq2[1])];

→ step 6: we → we1,BC inside the list:we1

we1:[]$
for i:1 thru length(we) do(
hh:[we[i][1],rhs(kp2[i]),we[i][2],we[i][3],we[i][4],we[i][5],rhs(kq2[i])],    
we1:endcons(hh,we1)
)$

we1$

step 7: graph → we1,BC inside the list:we1

grp1[i]:=sol_points(we1,first,i);

wxdraw2d(point_type=-1,
points_joined = true,color=red,grp1[second],color=blue,grp1[third],
color = green,grp1[fourth],color = orange ,grp1[fifth],color=black,
grp1[sixth],color =yellow,grp1[seventh], yrange = [-12,2],  
title="pde→differential solution(=w1 and BC),oc:method");

color:red:c[1](t),blue:c[2](t),green:c[3](t),orange:c[4](t),black:c[5](t),yellow:c6[t]
==================================================================
rem: reason negative values ( = see above ,wb1 → -wb1),reason:graph is upsite down.

wxdraw2d(point_type=-1,
points_joined = true,color=red,grp1[second],color=blue,grp1[third],
color = green,grp1[fourth],color = orange ,grp1[fifth],color=black,
grp1[sixth],color =yellow,grp1[seventh], yrange = [-0.7,0.7],  
title="pde→differential solution(=w1 and BC),oc:method");