Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am stacked with a partial differential equation as follows: This is the $\LaTeX$ expression for the equation

$$ \frac{\partial y}{\partial t}=\frac{\partial^2 y(x,t)}{\partial x^2} -\left(a-\frac{1}{\cosh ^2(x)}\right) y(x,t)-y^3(x,t) $$

The Mathematica code is here. This is the equation with with $a=0.2$ to be specific:

eq = D[y[x, t], t] == 
   D[y[x, t], {x, 2}] - 0.2 *y[x, t] + Cosh[x]^-2*y[x, t] - y[x, t]^3;

and here are the boundary and initial conditions:

bc1 = y[-10, t] == 0;
bc2 = y[10, t] == 0;
ic = y[x, 0] == 0.001*(Exp[-x^2] - Exp[-10^2]);

Here I solve and plot it:

s = NDSolve[{eq, bc1, bc2, ic}, y[x, t], {x, -10, 10}, {t, 0, 20}, 
   PrecisionGoal -> 2];

Plot[Evaluate[y[x, t] /. s /. t -> 20], {x, -10, 10}, 
 AxesLabel -> {Style["x", 16, Italic], Style["y(t=20)", 16, Italic]}]

The result looks like this:

enter image description here

and qualitatively it is what I expect.

It should be explained here that this equation has a fixed point. Depending upon the parameter $a$ it may be either $y=0$, or a solution like the one shown in the image above. I am only interested in the system behavior in the fixed point. Intermediate steps are of no interest. For this reason I plot the result at $t=20$, expecting that the solution is already close enough to the fixed point.

Let us check the conversion of the obtained solution. The simplest way to do that is to look at the evolution of $y(0,t)$. Here it is:

enter image description here

The conclusion is that what we see is, so far, non-converging. I checked up to $t=200$ with the output of the same type showing non-convergence. It seems to be so for any time duration.

However, I solved the same equation with Comsol, and there it nicely converges after only a small number of time steps (five or so). For technical reasons I would, nevertheless, prefer it solved with Mathematica, rather than with whatever else.

My question: do you see any workaround?

share|improve this question

2 Answers 2

up vote 5 down vote accepted

Actualy, for $a=1/5$ the solution converges. The convergence speed for relaxation process is sensitive to initial condition, and greately slows down near fixed point solution. Using default options for NDSolve the solution converges to attractor after $t\approx 50$

Plot[Evaluate[y[0, t] /. s], {t, 0, tmax}, 
 AxesLabel -> {Style["t", 16, Italic], Style["y(x=0,t)", 16, Italic]}]

enter image description here

but it causes NDsolve::eerr. PrecisionGoal -> 2 solves it but this is not the best solution. The better approach is to use spectral spatial discretization

mol = {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 16, 
    "DifferenceOrder" -> "Pseudospectral"}};

s = NDSolve[{eq, bc1, bc2, ic}, y, {x, -10, 10}, {t, 0, tmax=100}, 
 Method -> mol];

Plot[Evaluate[{y[x, 20], y[x, 30], y[x, 40], y[x, 50]} /. s], {x, -10, 10},
 PlotRange -> All, 
  AxesLabel -> {Style["x", 16, Italic], Style["y(x,t)", 16, Italic]}]

enter image description here

share|improve this answer
    
Thank you for the answer. You are right, it is the matter of maximal time. With the value of the parameter a=0.365, at which I checked up to tmax=200 without the converging shows convergence at tmax>400. May be it is, indeed, method-dependent, since in Comsol converging has been achieved at much smaller tmax. Please, explain what criteria you used to choose the method? Why this method is better (for it is not visible in the plot y(0,t), I compared them obtained by your method and by Automatic)? Thank you once more. –  Alexei Boulbitch Oct 29 '13 at 8:41
    
@AlexeiBoulbitch I've used spectral methods just because of their spectral convergence property for smooth functions, but You are free to choose FDA also (Automatic). –  mmal Nov 11 '13 at 16:20

I just adopted the suggestions prompted in the non-convergence error message.

mol2 = {"MethodOfLines",  "SpatialDiscretization" -> {"TensorProductGrid", 
 "MinPoints" -> 1000}};

and Mma find the solution.

Maybe this help.

share|improve this answer
    
Thank you, but the question arises, since on my machine (WinXP, Mma 9.0.1.0) Mma gives no error message and, therefore, no suggestion after my original code. It is, probably, just the place, where science transfers into art. –  Alexei Boulbitch Oct 29 '13 at 10:51

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.