Why is the PDE Toolbox solving my parabolic equation incorrectly?
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Catherine
el 11 de Nov. de 2013
Comentada: Catherine
el 18 de Nov. de 2013
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 comentarios
Marc
el 12 de Nov. de 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....
Respuesta aceptada
Bill Greene
el 12 de Nov. de 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 comentarios
Bill Greene
el 13 de Nov. de 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
el 13 de Nov. de 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
Más respuestas (1)
Bill Greene
el 18 de Nov. de 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
Ver también
Categorías
Más información sobre Eigenvalue Problems en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!