kill(all)$

maxima : appendix : fvm , Dirichlet BC applied (=left boundary),version 1.0
peter.vlasschaert@gmail.com,3/11/2022

ref :taylor,grid (CV)
------------------------
new:31/10/2022:maxima :  appendix : fvm : Grid generation for CV in 1D,version 1.0
new:20/10/2022:maxima :  appendix : fvm : first order ,finite difference (central,forward,backward)
new:20/10/2022:maxima :  appendix : fvm : second order ,finite difference (central),version 1.0

BC : Boundary Condition (used : left of the right of the interval)

load package,option
------------------------

load(lrats) $ ratprint:false$

representation : differential 
---------------------------------

s1:wx1=at(diff(Φ(x),x,1),x=w)$
s2:wx2=at(diff(Φ(x),x,2),x=w)$
s3:wx3=at(diff(Φ(x),x,3),x=w)$

s:[s1,s2,s3];

discritization : Taylor series 
----------------------------------

find : wx1, p=P (Grid : ref)

p1:Φ[p]=Φ[w]+Δx/2*wx1+1/2*(Δx/2)^2*wx2;
p2:Φ[E]=Φ[w]+3*Δx/2*wx1+1/2*(3*Δx/2)^2*wx2;

solve : doesn't work 
------------------------
solution : use eliminate "Maxima"

p3:eliminate([p1,p2],[wx2]);

p4:p3[1];

p4a:part(p4,0);

p4aa:part(p4,1,0);
p4bb:part(p4,1,1);
p4cc:part(p4,1,2);
p4dd:part(p4,1,3);

p5:p4dd = 0;

solve : wx1
--------------

p5a:solve(p5,wx1)[1];

check : [8-9+1 = 0]

p5a1:ratexpand(p5a);
p5a2:ratsimp((p5a));

p6b1:num(rhs(num(p5a)));

p6b2:coeff(p6b1,Φ[p]);

p6b3:coeff(p6b1,Φ[E]);

p7a:expand(p1*p6b2);
p7b:expand(p2.p6b3);

p7c1:p7a+p7b;

p7c2:solve(p7c1,wx1)[1];

check : wx1
--------------
? p5a,p7c2

is(rhs(p5a)=rhs(p7c2));

find : Error (wx3)
---------------------

p1a:Φ[p]=Φ[w]+Δx/2*wx1+1/2*(Δx/2)^2*wx2+1/6*(Δx/2)^3*wx3;
p2a:Φ[E]=Φ[w]+3*Δx/2*wx1+1/2*(3*Δx/2)^2*wx2+1/6*(3*Δx/2)^3*wx3;

p8a:expand(p1a*p6b2);
p8b:expand(p2a.p6b3);

p9a:p8a+p8b;

p9b:expand(solve(p9a,wx1))[1];

error : second order : O(Δx^2), thrid order differential.

p9b1:ratcoeff(rhs(p9b),wx3);

p9c:ratsimp(subst([wx3=0],p9b));

check :wx1 
-------------
p5a & p9c

is(rhs(p5a)=rhs(p9c));

final result:
--------------

p10:subst(s,p9b);

→ first part : error (Δx^2)

p10a:part(rhs(p10),1);

→ second part : wx1

p10b:length(rhs(p10));

p10b1:sum(part(rhs(p10),i),i,2,p10b);

p10b2:ratsimp(p10b1);

→ combine :

p11:lhs(p10) =p10a+p10b2;

wx1 : need inside p11

p11a:expand(lratsubst([reverse(p5a)],p11));

---------------------------------------------------------------------------------------