# Diffusion equation without time marching (but reversed)

Here is like in advection_global_reversed.m, except that we do this for the diffusion equation (which should be much more ill-behaved, and more sensitive to errors).

clear all; clf

% parameters
Nx=61; % number of gridpoints
Nt=60; % gridpoints in time
Lx=10; % domain length in x
Lt=4; % duration in time
xpos=5; % position of the final condition
mu=0.1;

% differentiation
[d.x,d.xx,d.wx,x]=dif1D('cheb',0,Lx,Nx);
[d.y,d.yy,d.wy,y]=dif1D('fd',0,Lt,Nt,5);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

% rename y to time derivative
D.t=D.y;
D.tt=D.yy;
t=y;
T=Y;

l.left=[l.cbl;l.left;l.ctl];
l.right=[l.cbr;l.right;l.ctr];
l.start=l.bot;
l.end=l.top;

% system matrices
A=D.t-mu*D.xx;
b=zeros(NN,1);


Here this is the “final guess”.

% final guess
f0=exp(-((X-xpos)/1).^2);


# Boundary conditions

The “left”, “right” and “end” boundaries should be like in the initial guess.

% boundary conditions
loc=[l.end; l.right; l.left];
C=I(loc,:);

b(loc)=C*f0(:);
A(loc,:)=C;


# Solve system

f=A\b;
f=reshape(f,Nt,Nx);

% show evolution of f
subplot(1,2,1);
surf(X,T,f); view(2)
xlabel('x'); ylabel('t');
title('time evolution of the string');


# March in time of the initial condition

Here this is some kind of validation, we use a different method to go back from the initial condition to the final one, and compare to the global solution.

% march in time the initial condition
I=eye(Nx);
E=I;
A=mu*d.xx;

% boundary conditions
loc=[1;Nx];
E(loc,:)=0;
A(loc,:)=I(loc,:);

% march in time matrix
dt=t(2)-t(1);
Mm=(E-A*dt/2);
M=Mm\(E+A*dt/2);

% initial condition
q=f(1,:)';

% marching loop
for ind=1:(Nt-1); q=M*q; end

% plotting
subplot(1,2,2);
plot(x,q,'b',X(l.end),f(l.end),'r.-');
axis([0,Lx,0,1.2])
legend('march in time','global'); title('diffusion reversed')
xlabel('x'); ylabel('final time');

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','diffusion_global_reversed.png');


The subplot on the left shows the evolution in time of the solution. You see that this is ill-behaved because the solution becomes very large and oscillatory.