Looking for a matlab/maple code for plotting the truncation error

Question

On page 18 on this text: http://www.dima.uniroma1.it/users/lsa_adn/MATERIALE/FDheat.pdf , the graph in figure 8 on this page, how would I write a suitable code in matlab or maple that will produce this graph? It's a graph of the different truncation errors of the different finite difference methods (on the last page of this text there are codes for the different finite difference schemes approximations to the exact solution.)

To my more immediate necessity, I have the PDE:

$$ u_t = u_{xx} +\sin(x+t)-\cos(x+t) , \\ u_x(0,t)=-\sin(t) , \\ u_x(1,t)=-\sin(1+t) , \\ u(x,0)=\cos(x)$$

and I want to find the truncation error of pdsolve's default method in maple and plot it in maple, I was told it should be a linear graph where the $Y$ axis is between $10^{-4}$ and $10^{-8}$. How to implement this in maple?

Answer 1

Just for creating the plot, here's the MATLAB code to generate a graph similar to the one you cited (I've mocked up the data)

Dx = 2.^-(2:8);

% These errors to be generated by your code: mocked up here
E_FTCS = Dx.^1.7;
E_BTCS = 0.9*Dx.^1.8;
E_CN = Dx.^1.75;

E_ideal = Dx.^2;

% Create the plot
loglog(Dx,E_FTCS,'ko',Dx,E_BTCS,'k*',...
       Dx,E_CN,'ks',Dx,E_ideal,'k--')

legend({'FTCS','BTCS','CN','ideal'},'Location','northwest')
xlabel('\Delta x')
ylabel('Error: ||u-u_e||_2')
title('\Delta t = 5.0e-6 (constant)')

Output graph with MATLAB version R2014b:

Answer 2

This can be done fairly easy in maple. For this I have choose the pde from your other question because that pde has an exact solution.

restart:with(plots):
pde:= diff(u(x, t), t) = diff(u(x, t), x, x)-sin(x+t)+cos(x+t);

First we will check whether your exact solution satisfy the pde or not. For this

ansol:=u(x, t)=cos(x+t);

0

The zero output means that the exact solution satisfies the pde. Now moving on to the numerical solution,

ic:={u(x,0)=cos(x)};
bcs:={ D[1](u)(0,t)=-sin(t),D[1](u)(1,t)=-sin(1+t)};

pds:= pdsolve({pde},ic union bcs,numeric,time=t,range=0..1,errorest=true,timestep=1/16,spacestep=1/16);

In pdsolve we have utilized the option errorest=true to compute error estimates.

The error estimators used, in both the visual error estimates and the error control, are simply local truncation error estimates for the PDE or PDE system.

First, lets find out the effect of the timestep and spacestep on the error,

pds:-settings(timestep=1/8,spacestep=1/8);
pds:-plot([u(x,t),[u(x,t)+err(u(x,t)),color=blue,linestyle=2],[u(x,t)-err(u(x,t)),color=green,linestyle=3]],t=5,axes=boxed);

Now we reduce the timestep and spacestep

pds:-settings(timestep=1/64,spacestep=1/64);
pds:-plot([u(x,t),[u(x,t)+err(u(x,t)),color=blue,linestyle=2],[u(x,t)-err(u(x,t)),color=green,linestyle=2]],t=5,axes=boxed);

Clearly, we can see from the above plot that the error in the pde solution is acceptable.

Now we will plot the error in the solution along the space

Plotting error in the solution with time,

Finally, visualizing the absolute error,

pds:-plot([[(abs(u(x,t)-(cos(x+t)))^2),color=red]],t=5,axes=boxed);

pds:-plot([[(abs(u(x,t)-(cos(x+t)))^2),color=red]],t=0..5,x=0.5,axes=boxed);