kill(all)$

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

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

load(lrats)$ ratprint:false$

ref : see website
=============

new: 17/11/2020:maxima :modified TRAM ODE-BVP problem ,version 1,part 1
download  (Maxima,code)       :maximamodifiedTRAMODEBVPversion1finalpart1
download  (Maxima,png)         :maximamodifiedTRAMODEBVPversion1finalpart1outp
new: 17/11/2020:maxima :appendix modified TRAM ODE-BVP problem ,version 1.
download  (Maxima,code)       :appendixmodifiedTRAMODEBVPversion1final
download  (Maxima,png)         :appendixmodifiedTRAMODEBVPversion1finaloutp

ODE : TRAM
===========
Pe : Peclet          ( = 6 )
Da : Damkohler  ( = 2)
Domain :  0 ≤ x ≤ 1

parameter 
---------------

Pe:6;
Da:2;

z1:1/Pe*'diff(y,x,2)-'diff(y,x)=Da*y^2 ;

ODE BVP : Boundary Value Problem. => BC

→ BC1:y[0]
=========

z1a:'diff(y,x)=Pe*(y-1);

BC1:y[N+2]
=========

z1b:'diff(y,x) = 0;

finite difference model :approximation model
=================
"isothermal  TRAM ,irreversible sec. order reaction ."

ODE : 
====

p[i]:=(y[i+1]-2*y[i]+y[i-1])/(Pe*dx^2)-(y[i+1]-y[i-1])/(2*dx) = Da*(y[i])^2;

p[1];

BC1:y[0] 
=======

q1:(y[2]-y[0])/(2*dx) = Pe*(y[1]-1);

q1a:solve(q1,y[0]);

BC1:y[N+2] = 'y[6],example'
=========

q2:(y[6]-y[4])/(2*dx)= 0;

q2a:solve(q2,y[6]);

BC 
===

q2b:flatten(cons(q1a,q2a));

number of equations
================

N:6 ; 

q2c:N-2; /*2+1,index start from zero ( =BC) ,from both Boundary Condition*/

q3:makelist(p[i+1],i,0,q2c);

rem:number of equations (q3) = q2c-0 +1

q4:subst(q2b,q3);

q4a:expand(q4);

q4b:makelist(y[i],i,1,q2c+1);

q4b1:apply(matrix,[q4b]);

vector: b
=======

q4b1a:transpose(q4b1);

matrix : A
========

q4b2:expand(coefmatrix(q4a,q4b));

q4b2a:expand(augcoefmatrix(q4a,q4b));

→  lhs : finite differences equations 

q4b2b:q4b2.q4b1a;

→  rhs : finite differences equations 

q5a:-1*col(q4b2a,6);

q5b:makelist(2*(y[i])^2,i,1,length(q4b));

q5ba:transpose(apply(matrix,[q5b]));

vector : c
=======

q6c:q5a+q5ba;

lhs = rhs
======

simp:false;

q7:q4b2.q4b1a = q6c;

simp:true;

q7a:q7;

length of interval : L
================
L : xb ≤ x ≤ xe , find : dx

xb:0;  /*xb  : xbegin from domain = L */
xe:1;   /*xe  : xbegin from domain = L */
L:xe-xb;
dx: L/(length(q4b)+1);

find eq : finite difference model
========================

q7b:ev(q7a,nouns); /*need : (evaluation 'ev') */

A*b = c , 'see above '
================

==================================
A    = matrix , 'square matrix => ' find  A^(-1)
b,c = column vector
==================================
I*b =A^(-1)*A*b = A^(-1)*c => b = A^(-1)*c
==================================

A = matrix 
========
use :q4b2

q8a:ev(q4b2,nouns);

A^(-1) = matrix 
===========

q8b:invert(q8a);

q8b1:float(q8b);

→ A.A^(-1) = unit matrix, dim(5) =unit(5)
================================

q8b2:q8a.q8b;

→ b = A^(-1)*c
============

start values : iteration process
=======================

q9:[1,0.8,0.6,0.4,0.2];

q9a:map("=",q4b,q9);

q9b:q4b1a=ev(lratsubst(q9a,q8b.q6c),nouns);

q9c:transpose(rhs(q9b));

convert : matrix ( = transpose(vector))  =>   list
----------------------------------------------------------------

→ new value (1)

q9d:flatten(args(q9c));

q9d1:map("=",q4b,q9d);

q9d2:q4b1a=ev(lratsubst(q9d1,q8b.q6c),nouns),float;

→ new value (2)

p10a:transpose(rhs(q9d2));
p10b:flatten(args(p10a));
p10c:map("=",q4b,p10b);
p10d:q4b1a=ev(lratsubst(p10c,q8b.q6c),nouns),float;

→ new value (3)

p11a:transpose(rhs(p10d));
p11b:flatten(args(p11a));
p11c:map("=",q4b,p11b);
p11d:q4b1a=ev(lratsubst(p11c,q8b.q6c),nouns),float;

→ new value (4)

p12a:transpose(rhs(p11d));
p12b:flatten(args(p12a));
p12c:map("=",q4b,p12b);
p12d:q4b1a=ev(lratsubst(p12c,q8b.q6c),nouns),float;

→ new value (5)

p13a:transpose(rhs(p12d));
p13b:flatten(args(p13a));
p13c:map("=",q4b,p13b);
p13d:q4b1a=ev(lratsubst(p13c,q8b.q6c),nouns),float;

→ new value (6)

p14a:transpose(rhs(p13d));
p14b:flatten(args(p14a));
p14c:map("=",q4b,p14b);
p14d:q4b1a=ev(lratsubst(p14c,q8b.q6c),nouns),float;

p15:ev(lratsubst(p14c,q7b),nouns),float;

example : How to use mnewton.
========================

load("mnewton")$

mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],
              [x1, x2], [5, 5]);

solve problem : finite difference
-------------------------------------------

1e) step : equation
--------------------------

p16a:q7b;

2e)  convert : matrix to list
------------------------------------

→ lhs : equations

p16a1:flatten(args(lhs(p16a)));

→ rhs : equations

p16a2:flatten(args(rhs(p16a)));

3e)  convert : matrix to list
------------------------------------

p16a3:map("-",p16a1,p16a2);

4e) list of variables
--------------------------

p16a4:listofvars(p16a3);

5e) initial values
--------------------------

p16a5:[0,0,0,0,0];

6e) solve
-------------

p16a6:mnewton(p16a3,
              p16a4, p16a5);

p16a7:flatten(p16a6);

p16a7a:rhs(p16a7[1]);

7e) check : use :p16a3
---------------------------------

p16a8:lratsubst(p16a7,p16a3),float;

8) make wxplot 2d
--------------------------

→ x values 

p16a9:makelist(dx*i,i,1,length(p16a5));

→ y values 

p16a10:makelist(rhs(p16a7[i]),i,1,length(p16a5));

→ (wx)plot2d, wx before 

l:[discrete, p16a9,p16a10];

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