2
$\begingroup$

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:

namely:

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)$

BCs:

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

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

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

and

Initial Value for s(x)=1;

$s \in[0,1]$;

My code looks like:

ClearAll["Global`*"];
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!

$\endgroup$
  • $\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$ – Henrik Schumacher Aug 22 '19 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$ – ABCDEMMM Aug 22 '19 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$ – ABCDEMMM Aug 22 '19 at 22:46
  • 1
    $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Aug 22 '19 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$ – Alexei Boulbitch Aug 24 '19 at 12:28
3
$\begingroup$

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

ClearAll["Global`*"];
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

$\endgroup$
  • $\begingroup$ why we use u[0, x] == u0 x/lr as boundary condition? $\endgroup$ – ABCDEMMM Aug 23 '19 at 15:22
  • $\begingroup$ @ABCDEMMM These are initial data, not boundary conditions. I took this to agree with the boundary conditions. $\endgroup$ – Alex Trounev Aug 23 '19 at 15:25
  • $\begingroup$ I donot define such initial value in the benchmark test. $\endgroup$ – ABCDEMMM Aug 23 '19 at 15:28
  • $\begingroup$ we donot use any body loading. just the displacement at the end of bar. $\endgroup$ – ABCDEMMM Aug 23 '19 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$ – Alex Trounev Aug 23 '19 at 15:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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