Based on the PDE model proposed by @Schumacher Solving a second order coupled PDE system, the one dimensional multi-field Problem I would like to solve such benchmark test:


u(x): displacement

s(x): diffusive parameter:

Strong forms:

$0=s^{2} E u^{\prime\prime}$

$0=s E (u^{\prime})^{2}-\mathcal{G}_{c}\left(2 \epsilon s^{\prime \prime}+\frac{1-s}{2 \epsilon}\right)$


$u(x=\pm L)=\pm u_{0}$

$s^{\prime}( \pm L)=0$

(Optimal:$s( \pm L)=1$)


Initial Value for s(x)=1;

$s \in[0,1]$;

My code looks like:

PDE1 = (s[x])^2 u''[x] == 0
PDE2 = 2 ϵ s''[x] + 0.5 (1 - s[x])/ϵ - 
   s[x] Ε (u'[x])^2/\[ScriptCapitalG] == 0

lr = 25;
ll = -25;
Ε = 1;
 \[ScriptCapitalG] = 1;
u0 = 10;
ϵ = 0.125;
ics = {u[ll] == -u0, u[lr] == u0, s'[ll] == 0, s'[lr] == 0, 
  s[ll] == 1, s[lr] == 1}
{uu,vv} = NDSolve[{PDE1, PDE2, ics}, {u, s}, {x, ll, lr}]

This code cannot run due to initial BCs problems in MMA.

DSolve also not possible in MMA:

DSolve::argm: DSolve called with 2 arguments; 3 or more arguments are expected.

Okay, then I have the following questions:

  1. can we solve this nonlinear problem in MMA11.3?

  2. how can we fix BCs in MMA, it should works? it seems not that difficult.

  3. can we use fdm or spectral method for NDsolve, or we have to use fem here?

  4. which numerical method is mostly recommonded for this problem (robust and efficient) for version 11.3?

Any help would be greatly appreciated

Many thanks in advance!

  • $\begingroup$ This is a system of two second order elliptic equations. You can only apply $2 \times 2$ boundary conditions: two for each elliptic equation. $\endgroup$ Aug 22, 2019 at 22:42
  • $\begingroup$ @HenrikSchumacher hey, thanks, in this case, I would like to get inhomogeneous solution, therefore I add two additional at the end of bar. $\endgroup$
    Aug 22, 2019 at 22:45
  • $\begingroup$ @HenrikSchumacher moreover, if I use the first four BCs, then I got the errors:NDSolve::berr: The scaled boundary value residual error of 4.948245340445189`*^77 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found. $\endgroup$
    Aug 22, 2019 at 22:46
  • 1
    $\begingroup$ You're welcome! $\endgroup$ Aug 22, 2019 at 22:52
  • 2
    $\begingroup$ The whole life I work with Ginzburg-Landau equation, but have never seen it in this form. Even more, I do not see nonlinearity in it. Could you please kindly comment, why do you call these equations Ginzburg-Landau? $\endgroup$ Aug 24, 2019 at 12:28

1 Answer 1


The method of the false transient works well to solve this problem. The solution quickly converges with an increase in the number of iterations t

PDE1 = (s[t, x])^2 D[u[t, x], x, x] == D[u[t, x], t];
PDE2 = 2 \[Epsilon] D[s[t, x], x, x] + 0.5 (1 - s[t, x])/\[Epsilon] - 
    s[t, x] \[CapitalEpsilon]1 (D[u[t, x], x])^2/\[ScriptCapitalG] == 
   D[s[t, x], t];

lr = 25;
ll = -25;
\[CapitalEpsilon]1 = 1;
\[ScriptCapitalG] = 1;
u0 = 10;
\[Epsilon] = 0.125;
bcs1 = {u[t, ll] == -u0, u[t, lr] == u0};
bcs2 = {s[t, ll] == 1, s[t, lr] == 1};
bcs3 = {Derivative[0, 1][s][t, ll] == 0, 
   Derivative[0, 1][s][t, lr] == 0};
ic = {u[0, x] == u0  x/lr, s[0, x] == 1};
{uu, vv} = 
 NDSolveValue[{PDE1, PDE2, bcs1, bcs2, ic}, {u, s}, {x, ll, lr}, {t, 
   0, 10}]

{Plot3D[uu[t, x], {x, ll, lr}, {t, 0, 10}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> Automatic], 
 Plot3D[vv[t, x], {x, ll, lr}, {t, 0, 10}, Mesh -> None, 
  ColorFunction -> Hue, PlotRange -> All, AxesLabel -> Automatic]}
{Plot[uu[10, x], {x, ll, lr}], 
 Plot[vv[10, x], {x, ll, lr}, PlotRange -> All]}

Figure 1 If you use the Newman condition, then the solution will be

{uu1, vv1} = 
 NDSolveValue[{PDE1, PDE2, bcs1, bcs3, ic}, {u, s}, {x, ll, lr}, {t, 
   0, 10}]

{Plot3D[uu1[t, x], {x, ll, lr}, {t, 0, 10}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> Automatic], 
 Plot3D[vv1[t, x], {x, ll, lr}, {t, 0, 10}, Mesh -> None, 
  ColorFunction -> Hue, PlotRange -> All, AxesLabel -> Automatic]}

Figure 2

  • $\begingroup$ why we use u[0, x] == u0 x/lr as boundary condition? $\endgroup$
    Aug 23, 2019 at 15:22
  • $\begingroup$ @ABCDEMMM These are initial data, not boundary conditions. I took this to agree with the boundary conditions. $\endgroup$ Aug 23, 2019 at 15:25
  • $\begingroup$ I donot define such initial value in the benchmark test. $\endgroup$
    Aug 23, 2019 at 15:28
  • $\begingroup$ we donot use any body loading. just the displacement at the end of bar. $\endgroup$
    Aug 23, 2019 at 15:29
  • 1
    $\begingroup$ @ABCDEMMM Firstly, you did not solve this problem. Therefore, it is impossible to call this a benchmark test. Secondly, when solving by the method that I proposed, it is necessary to set the initial conditions (in fact, this is the initial approximation to the solution). You can choose any initial conditions. I chose the ones that I liked. $\endgroup$ Aug 23, 2019 at 15:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.