kill(all)$

maxima : pde : general analysis advection pde :u't+c*u'x = 0,version 1.0,part 3
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
maxima : pde : general analysis advection pde :u't+c*u'x = 0,version 1.0,part 2

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

load(lrats)$

ref :

→ 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

von-Neumann Stability Analysis
-------------------------------------------

s1a:'diff(u,t) =sum(d[2*i+1]*'diff(u,x,2*i+1),i,0,v1)+sum(D[2*i]*'diff(u,x,2*i),i,1,v2);

part 1 :() ,  u =exp(qt)*exp(%i*k*x) ,from fourier series.

s1b:u =exp(q*t)*exp(%i*k*x);

→ first step: subst

s1c:subst([s1b],s1a);

v1:3;
v2:3;

s1c1:s1c;

s1c2:ev(s1c1);

s1c3:ev(%,nouns);

second step : solve for q

s2:solve(s1c3,q);

thrid step : q = a+%i*b
-----------------------------
a = real         part
b = complex part

see : Menu : Simplify → Complex simplification {realpart or imagpart }

s3:[q=a+%i*b];
s3a : a= realpart(rhs(s2[1]));
s3b : b= imagpart(rhs(s2[1]));

general : a ,equation:
---------------------------

s3a1:a='sum(D[2*j]*(-1)^j*k^(2*j),j,1,v11);

s3a2:a=subst([v11=v2],rhs(s3a1)),simpsum;

check :a

s3a3:is(rhs(s3a)=rhs(s3a2));

general :b ,equation:
---------------------------

s4a1:b='sum(d[2*j+1]*(-1)^j*k^(2*j+1),j,0,v11);

s4a2:b=subst([v11=v1],rhs(s4a1)),simpsum;

check :b

s4a3:is(rhs(s3b)=rhs(s4a2));

fourth step :stability

use : general ( use a,equation)

s5a:s3a1;

stability  : a ≤ 0,u=%e^(%i*k*x+q*t)),q = -a±%i*b
--------------------------------------------------------------

s5b:args(s5a);

→ j ≤ 1, j =integer.

s5c:args(s5b[2])[1]≤0;

for i:1 thru v2 do(
print("i=",i,"D[",2*i,"] =",subst([j=i,k=1],s5c))

)$

result :D[2] > 0,D[4] < 0,D[6]> 0...........

more general :
==========0

s6:'diff(u,t) = sum(r[i]*'diff(u,x,i),i,1,m);

use :s1b

s6a:s1b;

→ subst :

s6b:subst([s6a],s6);

s6b1:ev(s6b,nouns);

s6b2:subst([m=v1],s6b1);

s6b3:ev(s6b2,nouns);

solve :q

s6b4:solve(s6b3,q)[1];

find :r[i],i =1,2,3
----------------------

→ r[1]

s7a:'diff(lhs(s6b4),k,1) =diff(rhs(s6b4),k,1);

s7a1:at((s7a),[k=0]);

→ solve :

s7a2:solve(s7a1,r[1]);

→ r[2]

s7b::'diff(lhs(s6b4),k,2) =diff(rhs(s6b4),k,2);

s7b1:at((s7b),[k=0]);

s7b2:solve(s7b1,r[2]);

→ r[3]

s7c::'diff(lhs(s6b4),k,3) =diff(rhs(s6b4),k,3);

s7c1:at((s7c),[k=0]);

s7c2:solve(s7c1,r[3]);

s8:flatten([s7a2,s7b2,s7c2]);

result : use s6b4 

s9:subst(s8,s6b4);

find q : FTUS
-------------------

use : part 2( = s3b1,)

→ first step : given 

s10:(u[i]^(n+1)-u[i]^n)/Δt=(aa[1]*(u[i]^n-u[i-1]^n))/Δx;

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

use : s1b

s10a:s1b;

→ second step : discreet for s10a

general :u

s10b:u[i]^(n)=rhs(s10a);

→u[i+1]^(n):

s10b1:u[i+1]^(n) =subst([x=x+Δx],rhs(s10b));

→u[i-1]^(n):

s10b2:u[i-1]^(n) =subst([x=x-Δx],rhs(s10b));

s10bb:[s10b1,s10b2,s10b3];

→u[i]^(n+1):

s10b3:u[i]^(n+1)=subst([t=t+Δt],rhs(s10b));

s10bb:[s10b,s10b1,s10b2,s10b3];

→ third step : subst ⟶ s10, for u[i]^(n+1)

s11a:ratsimp(solve(s10,u[i]^n)[1]);

→ fourth step : subst : s10bb into s11a

s12a:expand(subst(s10bb,s11a));

s13a:expand(s12a/rhs(s10b));

→ fifth step : solve :s13a, for q

s14a:solve(s13a,q);

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