4
$\begingroup$

I am solving a system of three time-dependent partial differential equations on a circular domain, using the finite element method. However, in the outer part of the domain ($r>R$) the equations are different than in the inner part and the transition from the outer system of equations to the inner one has to be smooth.

To do this I defined the first set of equations ($EqUx1$, $EqUy1$, $EqT1$) and the second one ($EqUx2$, $EqUy2$, $EqT2$). And finally defined the actual system by multiplying the first set by the logistic sigmoid and the second one by $1-$ the logistic sigmoid:

dist[x_, y_] := (x^2 + y^2)^0.5; sig[x_, y_] := LogisticSigmoid[Slope*(dist[x, y] - Centre)]; EqUx = sig[x, y] * EqUx1 + (1 - sig[x, y]) * EqUx2; EqUx = sig[x, y] * EqUy1 + (1 - sig[x, y]) * EqUy2; EqT = sig[x, y] * EqT1 + (1 - sig[x, y]) * EqT2; Operator = {EqUx, EqUy, EqT}; bcs = DirichletCondition[{Ux[x, y, t] == 0., Uy[x, y, t] == 0., T[x, y, t] == Subscript[T, eq]}, True ]; ics = {Ux[x, y, 0] == 0, Uy[x, y, 0] == 0, T[x, y, 0] == Subscript[T, eq]}; pde = Operator == {0, 0, 0}; {Ux0, Uy0, T0} = NDSolveValue[{pde, bcs, ics}, {Ux, Uy, T}, {t, 0., Subscript[tt, last]}, {x, y} \[Element] mesh, DependentVariables -> {Ux, Uy, T}];

The sigmoid defined in the second line basically creates an inverted bump function in two dimensions (the Slope and Centre variables are used to manipulate the shape of the bump). The one dimensional equivalent would look like this:

sigR[r_] := LogisticSigmoid[Slope*(Abs[r] - Centre)];

This works nicely, however, I need a bump function that does not just converge to 1 and 0, but is exactly 0 before some radius value and exactly 1 after another. I have found some good examples (https://math.stackexchange.com/questions/101480/are-there-other-kinds-of-bump-functions-than-e-frac1x2-1), but they all produce errors while solving the equations numerically (finite element method) because Mathematica badly handles the points of transition to the exact value (0 or 1) even though they are infinitely differentiable (from the example in the link above).

So I am looking for a numerically-friendly bump function that has an exact value of 1 in the inner region and 0 in the outer region (so basically, a symmetrical smooth step function), can anyone help?

EDIT

These are the actual final equations that go into the solver (sorry for the greek letters):

System of equations

Improve this question
$\endgroup$
5
  • $\begingroup$ I think there can be different reasons for your errors. Maybe Mathematica does not handle your convex combination correctly. The most important terms are the highest derivative terms. If the highest derivative terms are the same in all equations (e.g. Laplacian) then I would isolate that somehow. But I do not know your equations. $\endgroup$ Commented Jul 22, 2022 at 15:31
  • $\begingroup$ Concerning your equations: For sig[x,y]==1 there are no derivatives left. Your interpolation therefore seems very singular. Try doing something similar in a simple example, say a simple ODE, may also fail. $\endgroup$ Commented Jul 22, 2022 at 15:55
  • $\begingroup$ A function you can try is (1 - Clip[(Abs[x] + h - 1)/(2 h), {0, 1}])^2 (1 + 2 Clip[(Abs[x] + h - 1)/(2 h), {0, 1}]), for a suitably tiny value of h. $\endgroup$ Commented Jul 22, 2022 at 15:59
  • 1
    $\begingroup$ I just wanted to stress again that sig itself may not the problem here, but the interpolation you use. Your inner equations are not differential equations but boundary conditions. You seem to be using the interpolation as some ad-hoc way to implement boundary conditions, but the problem may be mathematically ill-posed this way. If you have more of science background, note that your EqUx1 and your EqUx2 have different units, and that should give you pause. $\endgroup$ Commented Jul 23, 2022 at 16:30
  • $\begingroup$ Thank you @user293787, the point of using the sigmoid is to create Dirichlet boundary conditions that are gradual. The equations describe a specific kind of heat transfer and the boundary temperature cannot be abruptly fixed in this physical environment. The fact about different units shouldn't really be a problem since these equations are just different operators. But, you are right about the problem becoming singular (thanks for the comment), that is probably what causes the solver to fail. $\endgroup$ Commented Jul 24, 2022 at 22:20

2 Answers 2

Reset to default
3
$\begingroup$

If

f=Piecewise[{{Exp[-1/(1-x^2)],Abs[x]<1}},0];

is not ok, then perhaps you can try multiplying it by a Gaussian:

g=f*Exp[-4*x^2];

Also g is smooth everywhere, and identically zero for Àbs[x]>1, but it looks much nicer, which might translate into better numerical behavior. Of course, you should also try increasing the constant 4 in the exponent if that helps.

enter image description here

Code:

Plot[{f,g},{x,-2,2},PlotLegends->{"f","g"}]
Improve this answer
$\endgroup$
2
  • $\begingroup$ I am actually looking for a function that looks more like a box, for example: it gives exactly 1 for input -1 to 1, than there is a smooth transition so that in the regions before -2 and after 2 it gives exactly 0. $\endgroup$ Commented Jul 22, 2022 at 14:00
  • $\begingroup$ Thanks for the answer, it seems that as long as the definition of the function is piecewise the solver has a problem at the point of transition, even though it should be continuous. $\endgroup$ Commented Jul 22, 2022 at 14:57
1
$\begingroup$
  • Maybe use Cauchy function and the construction as in many differential manifold books.
cauchy[x_] = Which[x > 0, Exp[-1/x^2], x <= 0, 0]; g[x_] = cauchy[2 - x]/(cauchy[2 - x] + cauchy[x - 1]); bump[x_] = g[Sqrt[x*x]]; Plot[bump[x], {x, -3, 3}] Derivative[1][bump] /@ Subdivide[.5, 2.5, 20]

enter image description here

  • Mollifier should be another choice.
Clear[f, ϵ, ρ, φ]; ϵ = 1/10; f[x_] = Piecewise[{{1 + (1 + x + ϵ)/(1 - 2 ϵ), -2 + ϵ <= x <= -1 - ϵ}, {1, -1 - ϵ <= x <= 1 + ϵ}, {1 + (-1 + x - ϵ)/(-1 + 2 ϵ), 1 + ϵ <= x <= 2 - ϵ}}]; Plot[f[x], {x, -3, 3}] ρ = 1/NIntegrate[Exp[-1/(1 - x^2)], {x, -1, 1}]; φ[x_, ϵ_] = (ρ/ϵ)* Piecewise[{{Exp[-ϵ^2/(ϵ^2 - x^2)], -ϵ < x < ϵ}}]; Plot[NIntegrate[φ[t - x, .3]*f[x], {x, -2, 2}], {t, -3, 3}, PlotStyle -> Red, GridLines -> {{1, -1, 2, -2}, None}, GridLinesStyle -> Green]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for the answer, I tried implementing both methods and unfortunately, the calculations still fail. $\endgroup$ Commented Jul 23, 2022 at 14:17

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.