0
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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: enter image description here

$\endgroup$
  • $\begingroup$ Hi Steve, do you happen to know how do I solve my second question with the PDE that I posted in maple? Thanks. $\endgroup$ – Alan Jan 10 '16 at 17:52
  • $\begingroup$ I've not used maple in a while and I can't remember pdesolve: I think your solution is $u(x,t) = -\cos(x+t)$, so you somehow want to compare this with your output from pdesolve. Sorry I can't help more at the moment $\endgroup$ – Steve Jan 10 '16 at 18:16
  • $\begingroup$ I was wrong about the solution: it is if the signs are switched on your PDE. $\endgroup$ – Steve Jan 12 '16 at 11:49
0
$\begingroup$

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);

enter image description here

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);

enter image description here

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 enter image description here

Plotting error in the solution with time, enter image description here

Finally, visualizing the absolute error,

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

enter image description here

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.