# 1d coupled field PDE in MMA

how can we correctly solve 1d coupled field pde in MMA 12?

Description for the problem: the article the phase field fracture model in one dimension

I have proposed the following code (based on the solution from @Alex Trounev and @Cesareo), and found that:

Numerical solution from MMA Correct solution needs to look like: 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] (D[u[t, x], x])^2 == D[s[t, x], t];

lr = 25;
ll = -25;
u0 = 10;
\[Epsilon] = 0.125;
bcs1 = {u[t, ll] == -t u0 , u[t, lr] == t 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] == t u0 x/lr, s[0, x] == 1};
sols = NDSolve[{PDE1, PDE2, bcs1, bcs2, ic}, {u, s}, {x, ll, lr}, {t, 0, 10}]
utx = Evaluate[u[t, x] /. sols]
stx = Evaluate[s[t, x] /. sols]
Show[Plot3D[stx, {x, ll, lr}, {t, 0, 10}, Mesh -> None, ColorFunction -> Hue, AxesLabel -> Automatic, PlotRange -> All], PlotRange -> All]


PDE Model:

Coupled Eqs.:

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

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

BCs:

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

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

and

Initial Value for s(x)=1;

$$s \in[0,1]$$;

• Is that the equation you really want to solve? – Xminer Dec 26 '19 at 12:40
• @Xminer sure, see the updated infos for more details! – ABCDEMMM Dec 26 '19 at 14:16

Regarding the initial conditions

ic = {u[0, x] == t u0 x/lr, s[0, x] == 1};


there is an inconsistency because it is taken at $$t = 0$$.

Making

ic = {u[0, x] == 0, s[0, x] == 1};


it works fine.

sols = NDSolve[{PDE1, PDE2, bcs1, bcs2, ic}, {u, s}, {x, ll, lr}, {t, 0, 10}]
utx = Evaluate[u[t, x] /. sols]
stx = Evaluate[s[t, x] /. sols]
Show[Plot3D[utx, {x, ll, lr}, {t, 0, 10}, Mesh -> None, ColorFunction -> Hue, AxesLabel -> Automatic, PlotRange -> All],
Plot3D[stx, {x, ll, lr}, {t, 0, 10}, Mesh -> None, ColorFunction -> Hue, AxesLabel -> Automatic, PlotRange -> All], PlotRange -> All]
`
• I have changed the code, but it does not solve correctly, e.g.{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]} – ABCDEMMM Dec 26 '19 at 13:24
• @ABCDEMMM I included some additional commands to visualize the results. I used NDSolve instead NDSolveValue – Cesareo Dec 26 '19 at 13:45
• thanks you! but the solution from MMA is not correct, I have updated more details about this issue, see update infos! – ABCDEMMM Dec 26 '19 at 14:11