Why is the PDE Toolbox solving my parabolic equation incorrectly?

Question

Catherine 2013-11-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/105724-why-is-the-pde-toolbox-solving-my-parabolic-equation-incorrectly

评论： Catherine 2013-11-18

采纳的回答： Bill Greene

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Marc 2013-11-12

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 2013-11-12

在 MATLAB Online 中打开

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 2013-11-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/105724-why-is-the-pde-toolbox-solving-my-parabolic-equation-incorrectly#answer_114961

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

Bill Greene 2013-11-12

在 MATLAB Online 中打开

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

Bill Greene 2013-11-13

在 MATLAB Online 中打开

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 2013-11-13

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 2013-11-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/105724-why-is-the-pde-toolbox-solving-my-parabolic-equation-incorrectly#answer_115659

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Catherine 2013-11-18

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!

请先登录，再进行评论。

Why is the PDE Toolbox solving my parabolic equation incorrectly?

2 个评论
显示无隐藏无

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Why is the PDE Toolbox solving my parabolic equation incorrectly?

2 个评论 显示 无隐藏 无

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论