Discussion Forum

Generally my question relates to methods of coping with PDEs exhibiting bifurcations. I do it within the example of one rather simple equation.
To comprehensively formulate my question I need to start with a preamble. Please excuse me, if it is too long. I cannot explain my question otherwise.

1. Preamble.

1.1. Mathematical problem

I am solving a 2D stationary nonlinear equation exhibiting a bifurcation. The equation has the following form:

uxx+uyy=(A-sqrt(sqrt(x^2+y^2)+x)/sqrt(2*(x^2+y^2)))u+u^3 (1)

where u=u(x,y) is a dependent variable, A>0 is a parameter, and uxx+uyy is the Laplace operator applied to u. The problem should be solved within the infinite (x,y) plane with the Dirichlet boundary condition in infinity:
u(infinity)=0
Let me note that the lengthily expression sqrt(sqrt(x^2+y^2)+x)/sqrt(2*(x^2+y^2)) in (1) looks simple in the cylindrical coordinates, (r,fi), just as cos(fi/2)/sqrt(r).
At large, positive A values equation (1) has a stable trivial solution, u=0, while below some bifurcation value A=Ac>0 a nontrivial solution becomes stable, while the trivial one though still the solution, becomes unstable.

According to the bifurcation theory the value Ac may be found as the first eigenvalue of the auxiliary equation representing the linear version of the equation (1):

uxx+uyy=(Ac-sqrt(sqrt(x^2+y^2)+x)/sqrt(2*(x^2+y^2)))u (2)

I have an exact analytical solution for the equation (2) giving the bifurcation value: Ac=1/(2*2^1/3) approximately equal to 0.3969. I do not write it here, since it is rather long and I guess it will not help to answer the question that follows. I am ready, however, to give this analytical solution on demand, if anybody would find that it might help.

1.2. Numerical approach

I am trying to solve this equation with Comsol. What I need is the nontrivial solution u(x,y).
First I tried to solve Eq. (1) straightforwardly, as a stationary problem, and by fixing the boundary condition as u(R)=0, choosing R to be rather large, as, for example, R=100. This by far exceeds the feature size of the expected solution, which is about 10. I did it in the MathematicalEquations/General mode.
This did not work, since Comsol only returns the trivial solution, u=0, both at A>Ac and at A<Ac. I understand it in such a way that at the boundary one starts from the boundary condition u=0 which “sets the solver on the trivial solution”. So there is probably a mechanism to keep the solver at the trivial solution, even if it is unstable. I am not aware, if there is a way to switch this mechanism off.

I overcame this limitation introducing something like a pseudo-time stepping. That is, instead of solving the above stationary equation I am solving a pseudo-time dependent one:

du/dt - uxx-uyy=-(A-sqrt(sqrt(x^2+y^2)+x)/sqrt(2*(x^2+y^2)+eps))u-u^3 (3)

where t is the pseudo-time, and eps=0.00001 regularizes the equation in the point (0,0). Here I use the fact that the solutions of the stationary equation (1) I am looking for represent fixed points of the dynamic equation (3).
I solve (3) within a disk, with u(R)=0, up to the time T at which the solution becomes close to a stationary one. By trial I found that T=1000 is quite enough.
Please find the model file Bifurcation_Model.mph attached. In this file I removed all meshes and solutions. One may start with a Normal or Fine mesh size, provided the Adaptive Mesh Refinement is checked (which is the case in my model). The solution on my machine takes 2 munutes.

1.3. What such an approach gives and what it does not

Such a model, indeed, gives the trivial solution, u=0, at large A values, and the nontrivial solution emerges as soon as A becomes below some value AcNum. Here I introduce the bifurcation value AcNum, for the one obtained by simulation to distinguish it from the exact analytical value, Ac.
This can be checked by
(i) making the Parametric Sweep with respect to the parameter A from 0.37 to 0.43 with a step 0.005,

(ii) evaluating the Surface Maximum, Umax, in the Derived Values and

(iii) plotting the maximum value of the obtained solution, Umax, against A.

Please find the plot Umax=Umax(A) attached as the file Umax.png.
From the table corresponding to this plot I can derive the value AcNum as the A value at which Umax becomes greater than 0.0001.
The solution looks very much like what I expect.
Here, however, comes a problem.

2. Problem

The value of the parameter AcNum obtained this way is somewhat greater than the theoretical value. The difference, deltaAc=Ac-AcNum, not too great. deltaAc/Ac is about 5%. I could live with it, if only I understood
(i) where this difference comes from?
and
(ii) how it may be further decreased?
I need to know, how to get the AcNum close enough to the exact value Ac, to ensure that the solution is OK.
The difference, deltaAc, depends upon the choice of the mesh size. With the coarser mesh the AcNum decreases, which is expected Therefore, I tried
a) to decrease the mesh,
b) to use the complex meshes, such that the mesh is fine in the vicinity of the point (0,0) and is coarser far from it
c) to use the Adaptive Mesh Refinement
The outcome is counterintuitive: instead of coming closer to Ac from below, I find that AcNum goes to the value of about 0.422 that is higher than Ac. Any further mesh refinement does not help.
The best value deltaAc/Ac is about 5%.
My point is that no my action succeeded in decreasing this difference further, while some even increased it. I tried the following:
1) Varied R
2) Varied T
3) Changed relative and/or absolute tolerance
4) Used different solvers (not preconditioners though) and different time stepping mechanisms.
5) Defined a very fine mesh in the vicinity of the point x=0, y=0, with the Adaptive Mesh Refinement unchecked as well as with checked.
6) Solved the problem in cylindrical coordinates.

Nothing brings AcNum closer to Ac, while some of these actions make deltaAc/Ac even larger.
This is counterintuitive, especially the increase of deltaAc with the mesh size decrease. I would expect that the finer is the mesh; the better should be a solution.

I would like to stress that what I need is the solution u(x,y), rather than the bifurcation value AcNum. The latter I know from the analytical result. However, the problem described above makes me doubting the quality of the solution.

3. My questions
Do you see the way to get AcNum closer to Ac?

Where this difference comes from?

Is there any better way to numerically study equations exhibiting bifurcations?