restart;
with(PDEtools);
with(plots);
with(LinearAlgebra);
pde := diff(u(t,x),t) - d*diff(u(t,x),x,x);
IBC := [u(0,x)=f(x),u(t,-1)=0,u(t,+1)=0];
f := x -> exp(-8*x^2);
d := 1;
pds := pdsolve(pde,IBC,numeric,timestep=0.01,spacestep=0.01);
pds:-animate(t=0..1);
d := 'd':
pde;
N := 4;
stencil := diff(u(t,x),t) = add(alpha[i]*u(t+i*s,x),i=0..N)+Error;
Err := convert(convert(series(solve(stencil,Error),s,N+1),polynom),D);
vars := indets(Err) minus {s,t,x,seq(alpha[i],i=-N..N)};
Err := collect(Err,vars,'distributed');
sys := [coeffs(Err,vars)];
ans := solve(sys,{seq(alpha[i],i=0..N)});
subs(ans,stencil);
degree(convert(series(`leadterm`(convert(subs(ans,solve(stencil,Error)),D)),h,20),polynom),h);
GenerateStencil := proc(F,N,{orientation:=center,stepsize:=h,showorder:=true,showerror:=false})
local vars, f, ii, Degree, stencil, Error, unknowns, Indets, ans, Phi, r, n, phi;
Phi := convert(F,D);
vars := op(Phi);
n := PDEtools[difforder](Phi);
f := op(1,op(0,Phi));
if (nops([vars])<>1) then:
r := op(1,op(0,op(0,Phi)));
else:
r := 1;
fi:
phi := f(vars);
if (orientation=center) then:
if (type(N,odd)) then:
ii := [seq(i,i=-(N-1)/2..(N-1)/2)];
else:
ii := [seq(i,i=-(N-1)..(N-1),2)];
fi;
elif (orientation=left) then:
ii := [seq(i,i=-N+1..0)]:
elif (orientation=right) then:
ii := [seq(i,i=0..N-1)]:
fi;
stencil := add(a[ii[i]]*subsop(r=op(r,phi)+ii[i]*stepsize,phi),i=1..N);
Error := D[r$n](f)(vars) - stencil;
Error := convert(series(Error,stepsize,N),polynom);
unknowns := {seq(a[ii[i]],i=1..N)};
Indets := indets(Error) minus {vars} minus unknowns minus {stepsize};
Error := collect(Error,Indets,'distributed');
ans := solve({coeffs(Error,Indets)},unknowns);
if (ans=NULL) then:
print(`Failure: try increasing the number of points in the stencil`);
return NULL;
fi:
stencil := subs(ans,stencil);
Error := convert(series(`leadterm`(D[r$n](f)(vars) - stencil),stepsize,N+20),polynom);
Degree := degree(Error,stepsize);
if (showorder) then:
print(cat(`This stencil is of order `,Degree));
fi:
if (showerror) then:
print(cat(`This leading order term in the error is `,Error));
fi:
convert(D[r$n](f)(vars) = stencil,diff);
end proc:
GenerateStencil(diff(u(t,x),x$2),6,orientation=left,showerror=true);
centered_space := GenerateStencil(diff(u(t,x),x,x),3);
forward_time := GenerateStencil(diff(u(t,x),t),2,orientation=right,stepsize=s);
FTCS_stencil := subs(centered_space,forward_time,pde);
FTCS_stencil := expand(isolate(FTCS_stencil,u(t+s,x)));
Subs := [u(t+s,x)=u_future,u(t,x)=u_middle,u(t,x-h)=u_left,u(t,x+h)=u_right];
FTCS := unapply(rhs(subs(Subs,FTCS_stencil)),u_left,u_middle,u_right,d,s,h);
FTCS_solver := proc(f,tau,d,N,M)
local X, T:
X := j -> -1 + 2*j/(M+1);
T := i -> tau*i/N;
s := evalf(T(1)-T(0));
h := evalf(X(1)-X(0));
XX := Array(0..M+1,[seq(X(j),j=0..M+1)]);
u_past := Array(0..M+1,[seq(f(X(j)),j=0..M+1)],datatype=float);
u_past[0] := 0;
u_past[M+1] := 0;
#p[0] := plot(Matrix([[XX],[u_past]]^%T));
p[0] := plot(Matrix([[XX],[u_past]])^%T);
for i from 1 to N do;
u_future := Array(0..M+1,datatype=float);
u_future[0] := 0;
u_future[M+1] := 0;
for j from 1 to M do:
u_future[j] := FTCS(u_past[j-1],u_past[j],u_past[j+1],d,s,h)
od;
ArrayTools[Copy](u_future,u_past);
p[i] := plot(Matrix([[XX],[u_past]])^%T);
od;
display(convert(p,list),insequence=true);
end proc;
f := x -> x^2;
FTCS_solver(f,1,1,63,10);
pde;
centered_space;
MOL := subs(centered_space,pde);
MOL := subs(x=x[j],MOL);
Subs := [seq(u(t,x[j]+jj*h)=U[j+jj](t),jj=-1..1)];
MOL := subs(Subs,MOL);
U := 'U':
sys := 'sys':
M := 5;
U[0] := t -> 0;
U[M+1] := t -> 0;
for jj from 1 to M do:
sys[jj] := eval(subs(j=jj,MOL));
sys[jj] := isolate(sys[jj],diff(U[jj](t),t));
sys[jj] := rhs(sys[jj]) = lhs(sys[jj]);
print(sys[jj]);
od:
A := GenerateMatrix(convert(sys,list),[seq(U[j](t),j=1..M)])[1];
P[L] := 1 - s/2*A;
_A := (M,d,h) -> d/h^2*BandMatrix([1,-2,1],1,M);
_A(6,d,h);
CN := proc(f,tau,d,N,M)
local X, T:
X := j -> -1 + 2*j/(M+1);
T := i -> tau*i/N;
s := evalf(T(1)-T(0));
h := evalf(X(1)-X(0));
XX := Vector(M,[seq(X(j),j=1..M)],datatype=float);
u_past := map(x->f(x),XX);
p[0] := plot(Matrix([XX,u_past]));
P[L] := 1 - s/2*_A(M,d,h);
P[R] := 1 + s/2*_A(M,d,h);
for i from 1 to N do:
u_future := LinearSolve(P[L],P[R].u_past);
u_past := LinearAlgebra[Copy](u_future);
p[i] := plot(Matrix([XX,u_past]));
od:
display(convert(p,list),insequence=true):
end proc:
f := x -> x^2;
N := 100;
M := 10;
movie[1] := FTCS_solver(f,1,1,N,M):
movie[2] := CN(f,1,1,N,M):
display(Array([movie[1],movie[2]]));
TTdSMApJN1JUQUJMRV9TQVZFLzQzMDcxNDM2NTZYLCUpYW55dGhpbmdHNiI2IltnbCEiJSEhISM6IiYiJiwkKiYlImRHIiIiJSJoRyEiI0YsRigiIiFGLUYtRihGJ0YoRi1GLUYtRihGJ0YoRi1GLUYtRihGJ0YoRi1GLUYtRihGJ0YlTTdSMApJN1JUQUJMRV9TQVZFLzQzMTA0Njg0NDhYLCUpYW55dGhpbmdHNiI2IltnbCEiJSEhISM6IiYiJiwmIiIiRigqKCUic0dGKCUiZEdGKCUiaEchIiNGKCwkRikjISIiIiIjIiIhRjJGMkYuRidGLkYyRjJGMkYuRidGLkYyRjJGMkYuRidGLkYyRjJGMkYuRidGJQ==TTdSMApJN1JUQUJMRV9TQVZFLzQzMTA0Njg4MTZYLiUpYW55dGhpbmdHNiMmJSViYW5kRzYkIiIiRik2IltnbCEiJCEhISMzIiciJyIiIiIhLCQqJiUiZEdGKSUiaEchIiNGMEYtRi1GLEYtRi1GLEYtRi1GLEYtRi1GLEYtRi1GLEYrNiI=