kill(all)$

maxima : pde : general analysis advection pde :u't+c*u'x = 0,version 1.0,part 2
peter.vlasschaert@gmail.com,20/02/2023

ref :
maxima : pde : general analysis advection pde :u't+c*u'x = 0,version 1.0,part 1

load: packages.
--------------------

load(lrats)$

How to analyze FD scheme  , for linear PDEs
------------------------------------------------------------
FD    =  Finite Difference                    , numerical schemes
PDE =  Partial Differential Equations ,

→ 1e) advection & dispersion. ∂u/∂t = d[2n+1]*diff(u,x,2n+1),u=u(x,t)
      odd  : derivative
→ 2e) diffusion                        . ∂u/∂t = D[2n]*diff(u,x,2n)  ,u=u(x,t)
      even :derivatives

general pde equations  : n = 0,1,2,3......
-----------------------------------------------------

s1a:'diff(u,t) =d[2*n+1]*'diff(u,x,2*n+1);
s1b:'diff(u,t) =D[2*n]*'diff(u,x,2*n);

m = terms , rhs

s1:'diff(u,t) = sum(rhs(s1b)+rhs(s1a),n,0,m);

s1ca:subst([m=2,D[0]=0],s1),simpsum;

reason :D[0] = 0, only differential of x on rhs

general : m=5

s1cb :lhs(s1ca) = sum(aa[i]*'diff(u,x,i),i,1,5);

sc1b1:sort(args(rhs(s1cb)));

FD : scheme ,advection equation,see website : finite difference schemes
----------------------------------------------------------------------------------------------
i-1  → i → i+1 , distance(i,i-1) =distance(i+1,i) =Δx,distance(i-1,i+1) = 2*Δx

advection pde :

s2:lhs(s1cb) =sc1b1[1];

→ discreet , i = space , n = time
------------------------------------------
time  step :Δt,    n
space step:Δx,    i             see : above,distance

s2a1:u(x,t) =u[i]^n;

→ FTCS =( f)orward (t)ime & (c)entral differences (s)cheme

s2b1:[lhs(s2)=(u[i]^(n+1)-u[i]^n)/Δt, part(rhs(s2),2) =(u[i+1]^n-u[i-1]^n)/(2*Δx)];

use :s2

s2b2:subst(s2b1,s2);

method : unstable , unbounded

→ FTUS =( f)orward (t)ime & (U)pwind differences (s)cheme

first step : continue → discreet

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

→ flow direction 
                   <=     (up)wind for :( i+1 ) = i , index i , i =ref
i-1  → i → i+1

result :u[i+1]^n =u[i]^n , CS →UPS
------------------------------------------------
S  =scheme .
C   =central
UP =Upwind
------------------------------------------------

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

s3b1:subst([u[i+1]^n =u[i]^n,2*Δx =Δx],s2b2);

float(3!);

second step : taylor expansion
-----------------------------------------
f(x+Δx) =f(x)+Δx*f'(x)+(Δx)^2/2!*f'''(x)+...........  → u[i+1]
f(x-Δx) =f(x)-  Δx*f'(x)+(Δx)^2/2!*f'''(x)+...........   → u[i-1]

for i:1 thru 3 do(
print(i,"! = ",i!)
)$

→ find :u[i]^(n+1)  ,time taylor series

s3b2a:u[i]^(n+1) =u[i]^(n) +sum(Δt^i/i!*'diff(u[i]^n,t,i),i,1,m);

→ rhs :

s3part1:part(rhs(s3b2a),1);/*differential part*/

s3part2:part(rhs(s3b2a),2);

m1:3; /* finite taylor series : approximation */

s3b2b:subst([m=m1],s3b2a),simpsum;

s3b2b1:makelist(u[i],i,1,m1);
s3b2b2:subst([j=i],makelist(u[j],i,1,m1));

s3b2b3:map("=",s3b2b1,s3b2b2);

s3bt:subst(s3b2b3,s3b2b);

→ find :u[i-1]^(n),space taylor series

s4b2a:u[i-1]^(n) =u[i]^(n) +sum((-Δx)^i/i!*'diff(u[i]^n,x,i),i,1,m);

→ rhs

s4part1:part(rhs(s4b2a),1);/*differential part*/

s4part2:part(rhs(s4b2a),2);

s4b2b:subst([m=m1],s4b2a),simpsum;

s4bx:subst(s3b2b3,s4b2b);

thrid step :use taylor step 2 , into s3b1

s5:s3b1;

→ subst :{s3bt,s4bx}

s5a:ratsimp(subst([s3bt,s4bx],s5));

s5b:expand(s5a);

→ use : s2

s5c1:subst([part(lhs(s5b),m1) =0],lhs(s5b));

s5c2:s5b-s5c1;

s5c3:subst([Δx=0,Δt=0],s5c2);

s5c4:s5c2-rhs(s5c3);

→ replace : time derivatives → space derivatives.

=> first differentation :

s5t1:diff(lhs(s5c4),t)=aa[1]*'diff('diff(u[i]^(n),t),x);

s5t1a:'diff(u[i]^(n),t,2) =(aa[1])^2*'diff(u[i]^n,x,2);

→ general representation : wave equation (=s5t1a)

=> second differentation :

s5t1b:'diff(u[i]^(n),t,3) =(aa[1])^3*'diff(u[i]^n,x,3);

s5tt:[s5t1a,s5t1b];

→ time derivatives into s5c4

s6:subst(s5tt,s5c4);

s6a: lhs(s6) = w1*'diff(u[i]^n,x,2)+w2*'diff(u[i]^n,x,3);

s6w1:w1=factorsum(ratcoef(rhs(s6),'diff(u[i]^n,x,2)));

s6w2:w2=(ratcoef(rhs(s6),'diff(u[i]^n,x,3)));

s6w2a:lhs(s6w2) =aa[1]/6*v;

find : v

s6w2aa:solve(rhs(s6w2)=rhs(s6w2a),v);

s6w2:subst(s6w2aa,s6w2a);

use s6a

s7:subst([s6w1,s6w2],s6a);

we use . ∂u/∂t+c*∂u/∂x =0, book representation

s7a:subst([aa[1] =-aa[1]],s7);

s7b:subst([aa[1]=c],s7a);

s7c:truncate[errors] =rhs(s7b);

s7d:lhs(s7c) = diffusion(even :derivatives) +disperion(odd:derivatives)$

check : diffusion phenomena
I ) check coef : before even differentials in space.
II) question ;what positive diffusion ,solution boundd
                    what negativ diffusion  ,solution unbounded , blow up 

positive diffusion:
-----------------------

s8a1: Δx > c*Δt;
s8a2:1> c.Δt/Δx;

negative diffusion:
-----------------------

s9a1:Δx<c*Δt;
s9a2:1< c.Δt/Δx;

? what about : Δx = c*Δt, use :s7b
-------------------------------------------------------

same : modified finite difference equation : exact

p10a:[Δx=c*Δt];

p10b:subst(p10a,s7b);

p10c:solve(p10a,c);/*[m/s] =[m]/[s]*/

convergence :     c*Δt/Δx ≤ 1

p10d:subst(p10c,p10b);

what about :{Δx →0,Δt →0}
-------------------------------------
use:s7b

first → Δx = 0,second → Δt =0

p11a:limit(limit(lhs(s7b),Δx,0),Δt,0) ='limit('limit(rhs(s7b),Δx,0),Δt,0);

p11b:lhs(p11a) =ev(rhs(p11a),nouns);

when Δx→0,Δt→0 ,then result convergence to orginal eqaution,sse 
------------------------------------------------------------------------------------------
see above : truncate[error].

see :s2, "above "

→  exact equation

p12:s2;

p12a:lhs(p12)-rhs(p12) =0;

p12b:subst([aa[1]=-c],p12a);

→ when u = u[i]^n

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

rem : u[i+1]^n,taylor series in space .

use :s4b2a

p13a:s4b2a;

p13b:u[i+1]^n=subst([Δx=-Δx],rhs(p13a));

rem: s2b2, for x differential

p13c1:subst([m=m1],p13b),simpsum;

use :

p13c1:subst(s3b2b3,p13c1);

use :s4bx

p13c2:s4bx;

p13c3:p13c1-p13c2;

use :('diff(u[i]^n,x,3))=0

p13c4:subst([('diff(u[i]^n,x,3))=0],p13c3);

result : FD → Central Difference Scheme . first order differential with respect to space.
------------------------------------------------------------------------------------------------------------------

p13c5:reverse(p13c4/(2*Δx));/* see above :s2b1*/

explain :
-----------

n :time
i  : space

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