Why is the PDE Toolbox solving my parabolic equation incorrectly?

3 views (last 30 days)
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!
  2 Comments
Marc
Marc on 12 Nov 2013
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
Catherine on 12 Nov 2013
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!

Sign in to comment.

Accepted Answer

Bill Greene
Bill Greene on 12 Nov 2013
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
  5 Comments
Bill Greene
Bill Greene on 13 Nov 2013
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
Bill Greene on 13 Nov 2013
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

Sign in to comment.

More Answers (1)

Bill Greene
Bill Greene on 18 Nov 2013
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
  1 Comment
Catherine
Catherine on 18 Nov 2013
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!

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!