pdepe problem (numerical vs analytical solution)

So this is my pde:

The IC for 0<x<1:

The b.c for t>0.

u is bounded.

Using separation of variables by letting v = u-1, we then have the analytical solution:

The code in pdepe:

function [c,f,s] = pdex1pde(x,t,u,dudx)
c = 1;
f = dudx; 
s = 0;
end
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = 0; %ignored because m=2
ql = 0; %ignored because m = 2
pr = ur-1; 
qr = 0;
end
function u0 = u0(x)
u0 = 0; %initial condition
end

m = 2;
sol = pdepe(m,@pdex1pde,@u0,@pdex1bc,x,t);
u = sol(:,:,1);
x = linspace(0,1,50); %x = linspace(0,L,mesh points);
t = linspace(0,2,50);

Firstly, I do not see exponential decay or sinusoidal waves that should be expected from the analytical solution.

surf(x,t,u)
title('Numerical solution computed with 50 mesh points')
xlabel('Distance x')
ylabel('Time t')

Also, my analytical vs numerical solution does not match up.

Could someone tell me how to fix this?

plot(x,u(25,:),'o',x,1-(2*x.^(-1)*(pi)^(-1))'*exp(-pi^2*t(25))'*sin(pi*x))
title('Solution at t = 1')
legend('Numerical, 50 mesh points','Analytical','Location','South')
xlabel('Distance x')
ylabel('u(x,1)')

Answer 1

[Pat Gipper]

I don't think you were summing the terms of your series expansion correctly. The attached code carries out the expansion summation to 900 terms which seemed adequate. I also changed the source coefficient expression to what I thought is the correct equation. But regardless whether the source is your original, "s=0", or my answer, "s=(2/x) * dudx", there is still a mismatch between the series expansion and numerical solution. Oddly, I found that if you used "s=(1/x) * dudx" for the source there is very little error. Additionally, the results were no different whether a symmetry constant of m=1 or m=2 was selected.

In summary, I think the numerical solution is giving results that are making sense. We simply need to decide what source coefficient expression should be selected.