restart;
with(PDEtools);von Neumann stability analysisThe purpose of this worksheet is to introduce a few different stencils for the solution of the diffusion equation and to study their stability properties using the con Neumann stability analysis. We will use the centered_stencil and onesided_stencil procedures developed elsewhere:centered_stencil := proc(r,N,{direction := spatial})
local n, stencil, vars, beta_sol:
n := floor(N/2):
if (direction = spatial) then:
stencil := D[2$r](u)(t,x) - add(beta[i]*u(t,x+i*h),i=-n..n);
vars := [u(t,x),seq(D[2$i](u)(t,x),i=1..N-1)];
else:
stencil := D[1$r](u)(t,x) - add(beta[i]*u(t+i*h,x),i=-n..n);
vars := [u(t,x),seq(D[1$i](u)(t,x),i=1..N-1)];
fi:
beta_sol := solve([coeffs(collect(convert(series(stencil,h,N),polynom),vars,'distributed'),vars)]):
stencil := subs(beta_sol,stencil);
if (direction = spatial) then:
convert(stencil = convert(series(stencil,h,N+2),polynom),diff);
else:
subs(h=s,convert(stencil = convert(series(stencil,h,N+2),polynom),diff));
fi:
end proc:
onesided_stencil := proc(r,N,{direction := spatial})
local stencil, vars, beta_sol:
if (direction = spatial) then:
stencil := D[2$r](u)(t,x) - add(beta[i]*u(t,x+i*h),i=0..N-1);
vars := [u(t,x),seq(D[2$i](u)(t,x),i=1..N-1)];
else:
stencil := D[1$r](u)(t,x) - add(beta[i]*u(t+i*h,x),i=0..N-1);
vars := [u(t,x),seq(D[1$i](u)(t,x),i=1..N-1)];
fi:
beta_sol := solve([coeffs(collect(convert(series(stencil,h,N),polynom),vars,'distributed'),vars)]):
stencil := subs(beta_sol,stencil);
if (direction = spatial) then:
convert(stencil = convert(series(`leadterm`(stencil),h,N+1),polynom),diff);
else:
subs(h=s,convert(stencil = convert(series(`leadterm`(stencil),h,N+1),polynom),diff));
fi:
end proc:Using these procedures, we generate three specific stencils. The first two discretize a first order time derivative, while the last discretized a second order spatial derivative.forward_time := onesided_stencil(1,2,direction=temporal);
backward_time := subs(s=-s,forward_time);
centered_space := centered_stencil(2,3,direction=spatial);We will use these to develop a finite difference scheme for the diffusion problem:pde := diff(u(t,x),t)-d*diff(u(t,x),x,x);Our first scheme involves using the forward_time an centered space stencils. This gives the FTCS stencil we have previously considered.FTCS := subs(isolate(lhs(forward_time),diff(u(t,x),t)),isolate(lhs(centered_space),diff(u(t,x),x,x)),pde);If we instead use the backward_time stencil for the time derivative, we obtain the BTCS stencil.BTCS := subs(isolate(lhs(backward_time),diff(u(t,x),t)),isolate(lhs(centered_space),diff(u(t,x),x,x)),pde);
BTCS := subs(t=t+s,BTCS);Finally, taking the average of the FTCS and BTCS methods gives the Crank-Nicholson (CN) method. CN := (FTCS + BTCS)/2;It is interesting to compare these methods to the forward Euler, backward Euler, and trapezoidal method for solve the following ODE:ode := diff(y(t),t) = f(t,y(t));
forward_Euler := y(t+s) - y(t) - s*f(t,y(t));
backward_Euler := y(t+s) - y(t) - s*f(t+s,y(t+s));
trapezoidal := y(t+s) - y(t) - s*(f(t,y(t))+f(t+s,y(t+s)))/2;We see that the FTCS, BTCS, and CN methods are the PDE generalizations of the forward Euler, backward Euler, and trapezoidal methods. We can obtain the leading order error terms in each of the PDE methods by expanding each stencil in a power series in s and h and subtituting in the original PDE.pde_sub := isolate(pde,diff(u(t,x),t));
error_FTCS := convert(convert(dsubs(pde_sub,convert(series(series(FTCS,s,3),h,6),diff)),polynom),polynom);
error_BTCS := convert(convert(dsubs(pde_sub,convert(series(series(BTCS,s,3),h,6),diff)),polynom),polynom);
error_CN := convert(convert(dsubs(pde_sub,convert(series(series(CN,s,4),h,6),diff)),polynom),polynom);
Since the error in the FTCS stencil is O(s) + O(h^2), we call that method first order in time and second order in space, as is the BTCS method. But the CN stencil is second order in time and second order in space. Now, if we define the error E to be the difference between a numeric and true solutions for u(t,x), then it follows that E(t,x) obeys the following for each stencil: FTCS_E := subs(u=E,FTCS);
BTCS_E := subs(u=E,BTCS);
CN_E := subs(u=E,CN);The von Neumann stability analysis assumes the following ansatz for the error:E := (t,x) -> xi^(t/s)*exp(I*k*x);Putting this into each of E equations of motion and solving for the amplification factor xi, we obtainxi_FTCS := simplify(isolate(convert(simplify(expand(FTCS_E/E(t,x)),power),trig),xi));
xi_BTCS := simplify(isolate(convert(simplify(expand(BTCS_E/E(t,x)),power),trig),xi));
xi_CN := simplify(isolate(convert(simplify(expand(CN_E/E(t,x)),power),trig),xi));A particular method will be unstable if |xi|>0, as this implies that the error are growing exponentially in time. To determine under which circumstances this holds, we first simplify things by making the follwoing definitions:Subs := [s = -z*h^2/d/4,k=q/h];We arrange our solutions for xi into a list and make the above substitutionsxi_sols := simplify(subs(Subs,[xi_FTCS,xi_BTCS,xi_CN]));The extreme values of these are obtained when cos(q) = +1 or -1. Putting those values inextreme_xi_1 := subs(cos(q)=+1,xi_sols);
extreme_xi_2 := simplify(subs(cos(q)=-1,xi_sols));It is clear that extreme_xi_1 does not give us any useful information, but extreme_xi_2 does. We can now enforce the stability criterion |\134xi| < 1:stability_criteria_1 := map(x->subs(x,abs(xi)<1),extreme_xi_2);Notice these are exactly the same form as the stability criteria for the forward, backward and trapezoidal scheme for ODEs. In terms of the original parameters of the stencils, the stability criteria read:stability_criteria_2 := map(x->{rhs(x)<1,rhs(x)>-1},simplify(subs(isolate(Subs[1],z),extreme_xi_2)));These can be solved by MAPLE to give explicit conditions on s (assuming d is positive):stability_criteria_3 := map(x->solve(x,s),stability_criteria_2) assuming(s>0,d>0,h>0);The first element in the list says that the FTCS scheme is stable if s is less than h^2/2d, while the second and third elements say that the BTCS and CN methods are unconditionally stable.