Why is the PDE Toolbox solving my parabolic equation incorrectly?

Hello all,

I'm having a problem with getting the PDE Toolbox function 'parabolic' to return the correct solutions.

I'm working with the equation:

u_t = u_xx + u_yy +nu for any real c value

on the domain [-1, 1] x [-1, 1]

From the analytic solution, I know that with an initial condition of sin(pi*x)*sin(pi*y), the solution should be

u=e^(n-2pi^2)*t)*sin(pi*x)*sin(pi*y)

so if n < 2pi^2, then u should go to 0 over time. (e is raised to a negative exponent)

If n=2pi^2, then u should go to sin(pi*x)*sin(pi*y) (e is raised to 0)

and if n > 2pi^2, then u should go to positive infinity over time. (e is raised to a positive exponent)

I tried to simulate these results using PDE Toolbox and the parabolic function, but the results are far from what they're supposed to be as indicated by the analytic solution. Most of the solutions go to negative infinity, regardless of the n value. Can anyone tell me where I'm going wrong in my code/logic and why I'm getting wacky results?

I use the gui to define the boundaries and the mesh from x = -1:1 and y= -1:1, then I export decomposed geometry as p, e, and t.

Then, here's the MATLAB code I'm using to solve and plot:

tlist=0:.1:1; %Set tspan
c=1;
a=-(n); % I input actual values for n before I run the code.
f=0;
d=1;
u0='sin(pi*x).*sin(pi*y)';
u1=parabolic(u0,tlist,'squareb1',p,e,t,c,a,f,d);
pdemesh(p,e,t,u1(:,end))

The code runs and plots a solution, it's just not the correct solution according to the analytic solution.

*I'm using the raw code instead of the pde toolbox gui to solve simply because I prefer the raw code. I get the same results if I use the pde toolbox gui.

Any suggestions or comments are greatly appreciated!

[Marc]

I tried running your code and "p" is undefined. I get that you may have defined "p", "e" and "t" somewhere else, like in the GUI, but help me here try and help you....

[Catherine]

Yes, I've been defining my p, e, and t using the GUI. In the command window, I type pderect([-1 1 -1 1]), and the GUI window appears. From the top menu bar, I select mesh -> mesh mode to create a mesh, then mesh -> refine mesh to refine it. Then, I select mesh -> export mesh and p, e, and t are imported into my workspace.

Alternatively, this line of code has been working to define p, e, and t on [-1, 1] x [-1, 1]

[p,e,t]=initmesh('squareg'); % create a rectangular mesh 
[p,e,t]=refinemesh('squareg',p,e,t); %refine the mesh

Thank you so much for the help!

Answer 1

[Bill Greene]

Hi,

I have already replied to your post in the MATLAB newsgroup but I guess you didn't see it. I've included the same information below.

I tried your example in MATLAB R2013a with a = -.5 and the solution went to essentially zero at around .3 seconds.

I added these lines to the beginning of your code snippet to define the geometry and mesh.

gdm=[3 4 -1 1 1 -1 -1 -1 1 1]';
sf = 'S1';
nsm = ('S1')';
g = decsg(gdm, sf, nsm);
hmax = .2;
[p,e,t] = initmesh(g, 'hmax', hmax);

But that was my only significant change.

What version of MATLAB are you running?

Bill

[Catherine]

I'm using R2012a / Student Version 7.14

[Bill Greene]

Hi,

I tried this example on a 64-bit Windows version and a 64-bit Linux version of R2012a and I get the same result I got with the R2013a version. The solution in the interior is very close to zero at around .3 seconds.

I have attached my complete function below where I replaced my functions for mesh generation with the ones you posted. I have tried only a = -.5. Should I try a different value of a?

I forgot to ask the type of computer you are running on?

Bill

--------------------------------------------

function parabolic_problem_simple(  )
[p,e,t]=initmesh('squareg'); % create a rectangular mesh
[p,e,t]=refinemesh('squareg',p,e,t); %refine the mesh
tlist=0:.01:1; %Set tspan
c=1;
a=-(.5); % I input actual values for n before I run the code.
f=0;
d=1;
u0='sin(pi*x).*sin(pi*y)';
u1=parabolic(u0,tlist,'squareb1',p,e,t,c,a,f,d);
nid=getClosestNode( p, -.5, -.5 );
figure; plot(tlist, u1(nid,:));
figure; pdeplot(p,e,t,'xydata', u1(:,end))
figure; pdemesh(p,e,t,u1(:,end));
end

function nid=getClosestNode( p, x, y )
%getClosestNode Get node in point matrix p closest to input x and y
dist = (p(1,:) - x).^2 + (p(2,:) - y).^2;
[v, nid] = min(dist);
end

[Catherine]

I'm running it on a mac. I tried your code - it runs and n=0.5 goes to 0 as it's supposed to. When I get above 2*pi^2, though, the solution still goes off towards negative infinity instead of positive infinity. I noticed n=15 (should go to 0) shows some interesting behavior - it seems like the solution settles at 0 for a moment, then goes off in the negative direction. Is there anything else I should try?

Again, thank you so much!

[Bill Greene]

Hi,

I believe there are (at least) a couple of different things going on here.

First, if you use a finer mesh, e.g. something like this:

hmax = .05;
[p,e,t]=initmesh('squareg', 'hmax', hmax); % create a rectangular mesh

the case with a = -2*pi^2; will tend to infinity as time increases.

Until this morning, I hadn't looked closely at your analytical solutions and I still haven't looked closely enough. But are you sure about the factor of two in the exponential term (n-2*pi^2)? I'm wondering if that should simply be (n-pi^2)?

If I run the case a = -pi^2 it tends to zero-- until around five seconds and then begins to grow exponentially.

The numerical solution is obviously very unstable about the point u=e^0(...). It seems unlikely to me that it will be possible to compute this particular time-invariant solution.

Bill

[Bill Greene]

Hi,

I see my algebra was faulty-- so please ignore my comment about the factor of two. Sorry about that!

My comment about mesh density is still useful, however.

I'll continue to look at this.

Bill

Answer 2

[Bill Greene]

Hi,

I have looked further into this issue of a divergent solution when n is positive but less than 2*pi^2.

Not surprisingly, this behavior can be reproduced with the simpler 1D case. I get essentially the same results if I run the 1D case with the MATLAB pdepe function or with non-MathWorks PDE solvers.

The behavior can also be reproduced with a very simple numerical algorithm based on the central difference operator in the 1D spatial direction and a backward-Euler operator in time. It is not too difficult to show that the backward-Euler operator is not unconditionally stable when the n-coefficient is positive.

Numerical methods may exist for reliably solving this equation for positive n less than 2*pi^2 but I haven't been able to find them.

Bill

[Catherine]

The finer mesh seems to be helping with n values closer to 2*pi^2, but I've started to look into other numerical methods as well. Thank you so much for all of your help!