Message Boards Message Boards

0
|
8414 Views
|
8 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Please help to understand why NDSolve is not solving a nonlinear PDE in 2D?

Hello, I am solving the Belousov-Zhabotinsky system of PDE following the Gray-Scott scheme. When only the dynamical system is requested to be solved, NDSolve is giving the solution right away, however; when the diffusion term is added in 2D with Newman conditions at the boundaries (rectangular region), I tried different options and no one works. It seems to me that the numerical aspects of the problems (options for methods) are the key. What I expect to see? either concentric waves or spirals or maze-like patterns.

Attachments:
POSTED BY: David Quesada
8 Replies

No, I don't think that the non-linearity has something to do with the error.

And that is really mystic. I reduced the problem to one dimension AND removed the derivative of C2 in the boundary condition, but nevertheless the error message persists ???????

Look at the notebook.

Regards Hans

Attachments:
POSTED BY: Hans Dolhaine

Hi Hans: I found in the definition of the Laplacian in the second equation the index for the variable C it was C1 when it should be C2, cause the message of derivative of higher order. Second, I changed the way I wrote the Newmann boundary conditions at the frontier points. I fixed all of these and it is running now, however giving some errors of stiffness and that I should run with MinStepSize lower in order to capture the rapid variations. I did MinStepSize -> 0.005, takes a while for computing, it is not giving me yet what I am suppose to see, because the numerical error still persist, negative values for the concentration for example and exponential explosion of the numerical values, that does not make any sense. On the positive side, I have to say that the nonlinearity concern due to one of the functions is dissipating a bit.

POSTED BY: David Quesada

Ok. So the error message was helpful in finding the error in the pde's. Concerning the numerical problems I learned in former times, that sometimes a rescaling of the variables (concentration, time and lengths), so that constants have about the same order of magnitude, turns out to make things running smoothly. But I can't make a proposal having not looked intensely at your system.

Regards Hans

POSTED BY: Hans Dolhaine

miu = 0.0008; q = 0.55; eps1 = 0.04; kappa = 1; tmax = 100; (* time interval length *)

Lx = 1; (* x-side of the rectangle *)

Ly = 1; (* y-side of the rectangle *)

frontier = Rectangle[{-Lx, -Ly}, {Lx, Ly}]; C1ini = 1; C2ini = 1; freact[p, s] := p(1 - p) - ((p - miu)/(p + miu))q* s; (* reaction kinetics for eqn involving C1 *)

greact[p, s] := p - s; \ (* reaction kinetics for eqn involving C2 *)

eqn1 = D[C1[t, x, y], t] - D[C1[t, x, y], x, x] - D[C1[t, x, y], y, y] == freact[C1[t, x, y], C2[t, x, y]]/ eps1; (* diffusion-reaction eqn 1 *)

eqn2 = D[C2[t, x, y], t] - kappa(D[C1[t, x, y], x, x] + D[C1[t, x, y], y, y]) == greact[C1[t, x, y], C2[t, x, y]]; ( diffusion-reaction eqn 2 *)

InitialConditions = {C1[0, x, y] == C1ini, C2[0, x, y] == C2ini}; (* initial homogeneous concentration in the petri-dish *)

BoundaryConditions = {(D[C1[t, x, y], x] /. x -> -Lx) == (D[C2[t, x, y], x] /. x -> -Lx) == (D[C1[t, x, y], x] /. x -> Lx) == (D[C2[t, x, y], x] /. x -> Lx) == 0, (D[C1[t, x, y], y] /. y -> -Ly) == (D[C2[t, x, y], y] /. y -> -Ly) == (D[C1[t, x, y], y] /. y -> Ly) == (D[C2[t, x, y], y] /. y -> Ly) == 0}; \ (absence of flux through the boundaries ) solution = NDSolve[Join[{eqn1, eqn2}, InitialConditions, BoundaryConditions], {C1, C2}, {t, 0, tmax}, {x, -Lx, Lx}, {y, -Ly, Ly}]

Would the fact that freact[C1,C2] is a nonlinear function, create the problem for NDSolve? This is an example of partially nonlinear diffusion equation.

POSTED BY: David Quesada

Hmmmm - I boiled it a bit down and get an error message that boundary conditions and equations are not compatible. It says that the order of the derivatives in the boundary conditions must be smaller than in the equations. This is in my opinion the case, but....

Look at the notebook.

Regards Hans

Attachments:
POSTED BY: Hans Dolhaine

Hans, thank you once again. The derivatives by coordinates "x" and "y" are first order while the Laplacian is second order. I am just writing Newmann boundary conditions different from NewmannValue. I wrote then explicitly. Because the reaction term mixes C1 and C2, I am not sure if NewmannValue would be appropriate in that case.


Hans, thank you for your time. Here is the block dedicated only to the reaction diffusion solution. Comments were added to the end of each line to be understood. The workflow is as follows: Definition of the region (rectangular region) Definition of the reaction functions for each specie C1 and C2 Definition of the Diffusion operator (Dt C1 - Laplacian[C1]==f[C1,C2]) eqn1=Derivative by time of C1 - Laplacian2D[C1]== reaction function f[C1,C2] eqn2=Deivative by time of C2 - Laplacian2D[C2]==reaction function g[C1,C2] NDSolve is called for solution.

POSTED BY: David Quesada

Hello - the content of your notebook looks quite complicated. Could you give an abbreviated form which points out your problem?

Regards Hans

POSTED BY: Hans Dolhaine

Hans, thank you for your time. Here is the block dedicated only to the reaction diffusion solution. Comments were added to the end of each line to be understood. The workflow is as follows: Definition of the region (rectangular region) Definition of the reaction functions for each specie C1 and C2 Definition of the Diffusion operator (Dt C1 - Laplacian[C1]==f[C1,C2]) eqn1=Derivative by time of C1 - Laplacian2D[C1]== reaction function f[C1,C2] eqn2=Deivative by time of C2 - Laplacian2D[C2]==reaction function g[C1,C2] NDSolve is called for solution.

Attachments:
POSTED BY: David Quesada
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract