kill(all)$

maxima : fvm :appendix: linearized exponential scheme (=hybrid scheme ,a[W]/D[w]) part 2,version 1.0
peter.vlasschaert@gmail.com,16/12/2022

load package:
-----------------

load(lrats)$

→ see : website

ref 1:
new:15/12/2022:maxima : appendix : fvm : exponential scheme "1D -convection-diffusion,u > 0 ",version 1.0
use : {p10,p10a1,p10b1,p10c1},{p6b}
ref 2:
new:12/11/2022:maxima : :convection -diffusion relation peclet cell number (max :limit),version 1.0
find : limit =2
ref 3:
new:15/12/2022:maxima : fvm :appendix: linearized exponential scheme (=hybrid scheme ,a[E]/D[e]) part 2,version 1.0

hybrid scheme (=linerized exponential scheme).
====================================

use : general equation "exponential scheme"
----------------------------------------------------------

step 1:→ intro ref 1:general equation

p1a:a[p]*Φ(p) = a[W]*Φ(W) +a[E]*Φ(E);
p1b:a[p]=F[w]/(%e^Pe(Δx(w))-1)+F[e]·%e^Pe(Δx(e))/(%e^Pe(Δx(e))-1);
p1c:a[W]=F[w]·%e^Pe(Δx(w))/(%e^Pe(Δx(w))-1);
p1d:a[E]=F[e]/(%e^Pe(Δx(e))-1);

step 2:→intro ref 1: definition Pe(.), Peclet Cell Number

p1e:Pe(m1) = F[m]/D[m];

How get : Hybrid scheme (from :exponential scheme)
----------------------------------------------------------------------

make function :
===========

step 1: use "p1d"
---------------------

p2:p1c;

step 2: make p2 , "non-dimensional .a[W]"
----------------------------------------------
use : p1e

a) interface @ CS =[w]

p2a:[m1=Δx(w),m=w];

→ into p1e

p2a1:subst([p2a],p1e);

p2a2:[reverse(p2a1)];

→ divide p2 by D[e]

p2a3:p2/D[w];

→ use : p2a2 ,result : non-dimensional form

p2a4:lratsubst(p2a2,p2a3);

step 3: make function rhs(p2a4)
-----------------------------------------

→ set x = Pe(Δx(e)) =x

p3a1:subst([Pe(Δx(w)) =x],rhs(p2a4));

f(x) :=''p3a1;

rem : x-axis : Pe(Δx(e)) ,y-axis :a[E]/D[e] =f(x)

analysis function : f(x)
=============

step 1: make plot ,range : -5≤x≤ 5 for f(x)
-----------------------------------------------------

Menu :Plot→ Plot 2d, use : f(x), x =Pe(.)

p4:wxplot2d([f(x)], [x,-5,5],[y,-1,5], [legend, "f(x)t"],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

infinity limits
-------------------

p4a:'limit(f(x),x,inf) =limit(f(x),x,inf);

p4b:'limit(f(x),x,-inf) =limit(f(x),x,-inf);

rem : use L'Hopital rule : ∞/∞,0/0 = >  dif(num(f(x))/denom(f(x)))

limit @ x=0
-----------------

p4c:'limit(f(x),x,0) =limit(f(x),x,0);

linear:approximation around (x=0)
----------------------------------------------
linear : y =a*x+b constant ={a,b}

1e) method (Maxima : use taylor expansion)
---------------

Menu : Calculus :Get series:
Expression : f(x)
Variable     : x
Point          : 0
depth         : 8

taylor(, x, 0, 8);

p4d:taylor(f(x), x, 0, 8);/*infinite number of terms */

? How taylor be truncated some : depth  "first order"

p4d1:taylor(f(x), x, 0, 1);

p4d2:reverse(ratexpand(p4d1));

p4d3:f1(x):=''p4d2;

p4:wxplot2d([f(x),f1(x)], [x,-5,5],[y,-1,5], [legend, "f(x)t","f1(x)"],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

2e) method : differentation method (f(x)) with respect, line
-----------------------------------------------------------------------------

→ first step : differentation f(x) @ x=0

p5a1:diff(f(x),x);

rem : error when ,subst([x=0],%);

p5a2:'at('diff(f(x),x),x=0)='at(''diff(f(x),x),x=0);

problem : normal evaluation,use limit approach 

p5a3:lhs(p5a2) =limit(p5a1,x,0);

→ second step : use definition "tangent line (x0,y0)"

p5a4:y-y0=lhs(p5a3)*(x-x0);

x0 =0
y0 =use p4c

p5a5:subst([x0=0,y0=rhs(p4c),p5a3],p5a4);

→ third step : solve for y

p5a6:expand(solve(p5a5,y)[1]);

check :p4d2 ? p5a6
-------------------------

p6:is(p4d2=rhs(p5a6));

a) find intersection : y=0  & y=1+x/2
-----------------------------------------------

p7:subst([y=0],p5a6);
p7a:solve(f1(x)=0,x);

b) find value : blue curve  (= f(x)) when x=2, y1-value
----------------------------------------------------------------------

p8:y1 =subst(p7a,f(x));

p8a:ev(p8,numer);

c) what about large positive values of x for "f(x)"
----------------------------------------------------------------

p9:'f(1000) =f(1000);

p9a:lhs(p9) =ev(rhs(p9),number);

approximate 

p9a1:'f(1000)=(1000*%e^1000)/(%e^1000);

follows  x = 1000, y=1000

p9b:y = x;

p4:wxplot2d([f(x),f1(x),x], [x,-5,5],[y,-1,5], [legend, "f(x)t","f1(x)","x"],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

d) intersection of red & green line
---------------------------------------------

p10a:p5a6; /*red line*/
p10b:p9b; /*green line*/

→ intersection:

p10c:rhs(p10a) =rhs(p9b);

→ x-value : intersection

p10d:solve(p10c,x);

→ y-value : intersection

p10e:subst(p10d,p10a);

p4:wxplot2d([f(x),f1(x),x,0], [x,-5,5],[y,-1,5], [legend, "f(x)t","f1(x)","x","0"],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

linear profile : Pe(.)
--------------------------
------------------------------------------------------------------
lowest range : Pe =-2
→ d) intersection of red & green line, "see : above"
highest range : Pe =2
→ a) find intersection : y=0  & y=1+x/2,"see : above"
-------------------------------------------------------------------
see : ref 2
-------------------------------------------------------------------
Pe  > 2:
→ c) what about large positive values of x for "f(x)"
Pe < - 2
→ limit(f(x),x,inf) =0, "see:above,p4a"

Hybrid Scheme.
============

result :a[E]/D[e] 

a[W]/D[w] =0           if Pe <-2
a[W]/D[w] =1+Pe/2  if -2≤Pe≤2, "see above : method 1:p4d3,method 2:p5a6,x=Pe",CDS
a[W]/D[w] =Pe         if Pe>  2

result : combination upwind scheme, CDS scheme
-------------------------------------------------------------------
rem : CDS = Central Difference Scheme.

combination : combine with max function

p11:a[W]/D[w] =max(0,1+Pe/2,Pe);

rem : combination of straight lines , for Pe interval

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

point of intersection : y =1 ,both curve ( mirror curves,around y=1)
------------------------------------------------------------------------
----------------------------------------------------------------------------------------
blue curve  : part 1 → a[E]/D[e]    = max(Pe,1-Pe/2+1,0) 
blue curve  : part 2 → a[W]/D[w] = max(0,1+Pe/2,   Pe )
-----------------------------------------------------------------------------------------