I'm working in Mathematica to plot a 2-d catenary in Euclidean space according to this Math Stack Exchange Link, specified with a square Dirichlet boundary condition.
Here's what I've coded using Mathematica, followed by an error I receive. My goal was to set the surface integral to 1.5 because it's gotta be at least 1 if it's a sagging surface pinned at edges on a $1\times 1$ square.
NDSolve[
{ (-(D[u[x,y],x,x]+D[u[x,y],y,y])*(u[x,y]+λ)) == 1
, Integrate[
( 1+D[u[x,y],x]+D[u[x,y],y] )^(1/2) ∈ RectangularPolygon[4],x,y]==1.5
, DirichletCondition[u[x, y] == 0, True]
}
, u
, {x, y} ∈ RegularPolygon[4]
]
NDSolve: There are fewer dependent variables, {u[x,y]}, than equations, so the system is overdetermined
I think this is conditionally stable, so I've worked to provide a fair approximation of a stable error. If I omit the surface integral line, I'm told that NDSolve doesn't deal with non-linear coefficients. I was told elsewhere that this nonlinearity error is a hint that more detail is needed.
OK, I took a stab at improving the translation. I've got a coefficient array error that NDSolve can't parse, but it's closer to the problem statement. Surface integral not included yet.
NDSolve[
{
Sqrt[1+Norm[Grad[u[x,y],{x,y}]]^2 ]*Div[Grad[u[x,y],{x,y}]/Sqrt[1+Norm[Grad[u[x,y],{x,y}]]^2 ] ,{x,y}]==1/(u[x,y]+1),
DirichletCondition[u[x, y] == 0, True]
},
u,
{x, y} ∈ RegularPolygon[4]
]
The trouble with this problem is that a physical analog is basically impossible. Any two-dimensional approximation of a heavy surface, such as a sheet of chain male, introduces anisotropy and also costs a lot. A sheet of chamois was used to get a first stab, but, as shown, it seems to have something different going on in the corners.
Work done by @xzczd below shows a shape more like the Laplacian over a similar domain, possessing the similar, albeit less-pronounced, re-entrant geometry along a geodesic from the peak to the vertex (see here; Google image result).
What I'll probably need to do is generate a .CSV file with this and then hand-prototype a form from which I can cast a concrete part. CNC foam is a nice option for a daydream, but it costs about \$1500 x 2 where I'm at in the US Midwest. Rather than shell out the cash on something that might melt at cone 10 (pottery term), I may be better off measuring a grid of sticks, draping plaster gauze cloth over the peaks and smoothing subsequent dollops of hydrocal thereover with an eye on a printout to capture detail. I feign no skill at anticipating a successful method here, as I've not worked on something this big before.
It could come in handy that this surface looks like the native "dome" feature offered in SolidWorks, but, without sufficient documentation regarding the algorithm therein used, and also with there already being discrepancies between the Laplacian already looking 10% off from the @xzczd, I'm likely to err on the side of extra elbow grease and precision. If we've got the right math here, why use anything else?
A note on the picture above; chamois cloth (used wet as shown) stretches as it's comprised of, basically, an epidermal sponge, so there's probably some heat energy (and a smidge of elastic energy) absorbing the gravitational potential that deprives the final shape of any proximity to a more truthful catenary surface.
Dope work, all. Thanks to @xzczd.
@xzczd - This is a 38% scale prototype I'm making. If I can find affordable CNC service I'll use it, otherwise I'll make a pantograph and scale this up using a drill holding a ball mill following a ball stylus that runs over this form like a cam/template.
The ribs shown were generated using your Mathematica code modified a bit to output a $10 \times 10$ grid of points. I'm working with half of the shape at first (which is what you see). The next step will be to add the remaining webbing-like shapes that I've made from card stock.
Thank you - I think this is the first time in the history I've learned that I've seen this shape be made. I hope it works.
Been seeing similar shapes on domes for a few weeks now so doubt the immediate above.
Progress pics below. One shows a gauze cloth attempt that failed, but the next shows how I got to make the slip cast mold from solid plaster (similar in concept to the gauze cloth but without the cloth—two thin shells were made by pouring nearly-cured plaster slurry over the yellow papier mache and polyurethane master, then those two shells were joined with another thick overlay of plaster).
Current state, leather hard. Slip cast the main body and added hand built clay handles. Drying slowly now then firing and then testing. God help us all with the picture upload.
It’s probably 50-80% accurate. One big issue will be that the analytical surface is what I started with for the mold—the final shell thickness isn’t built in. If I get another chance I’ll do it, but for this moment I just want to know if I can cast something this large and then transport it to the field in one piece. If not, I may have to use castable refractory or try another approach.
Answer
I won't say it's easy to solve.
"I'm told that NDSolve doesn't deal with non-linear coefficients. " Yes, currently NDSolve cannot deal with nonlinear stationary PDE, so let's turn to finite difference method (FDM) to discretize $$\sqrt{1+|\nabla u|^2}\ \nabla\cdot{}\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}=\frac{1}{u+\lambda}$$. I'll use pdetoae for the generation of difference equations.
To take the constraint $$\int_{\Omega}\sqrt{1+|\nabla u|^2}dx=A$$ into consideration, I'll use trapezoidal rule to discretize it to a algebraic equation and solve it together with those difference equations.
points = 24;
{bL, bR} = domain = {0, 1};
grid = Array[# &, points, domain];
grad = Grad[#, {x, y}] &;
div = Div[#, {x, y}] &;
A = 3/2;
With[{u = u[x, y]}, sqrt = Sqrt[1 + grad[u].grad[u]];
eq = sqrt div[grad[u]/sqrt] (u + λ) == 1;
{bc@x, bc@y} = Table[u == 0 /. var -> b, {var, {x, y}}, {b, domain}]];
sqrtfunc = Function[{x, y}, #] &@sqrt;
difforder = 2;
(* Definition of pdetoae isn't included in this post,
please find it in the link above. *)
ptoafunc = pdetoae[u[x, y], {grid, grid}, difforder];
del = #[[2 ;; -2]] &;
ae = del /@ del@ptoafunc[eq];
aebc@x = del /@ ptoafunc[bc@x];
aebc@y = ptoafunc[bc@y];
rule = Flatten /@ (Outer[sqrtfunc, grid, grid] -> ptoafunc[sqrt]) // Thread;
trap[value_, grid_: grid] := 1/2 Differences[grid] (Most@value + Rest@value) // Total;
int = trap[Function[y, trap[Function[x, sqrtfunc[x, y]] /@ grid]] /@ grid] /. rule;
solfuncguess[x_, y_] := 0; λguess = 1/3;
solrule = FindRoot[{ae, aebc /@ {x, y}, int == A},
Append[Flatten[#, 1] &@
Table[{u[x, y], solfuncguess[x, y]}, {x, grid}, {y,
grid}], {λ, λguess}], MaxIterations -> 1000]; // AbsoluteTiming
(* {35.472657, Null} *)
solfunc = ListInterpolation[
Partition[Most[solrule][[All, -1]], points], {domain, domain}];
solrule // Last
(* λ -> 0.555491 *)
NIntegrate[sqrt /. u -> solfunc, {x, bL, bR}, {y, bL, bR}]
(* 1.51602 *)
Plot3D[{solfunc[x, y]}, {x, 0, 1}, {y, 0, 1}]
The achieved $A$ value isn't perfect, to achieve a better $A$, better initial guess is probably needed, which seems not to be that easy to find for this problem.
Guys read my relevant posts before (Well, will there be any?) may feel surprised that I'm not using Gauss-Legendre quadrature formula for numeric integration this time, that's because, with Gauss-Legendre quadrature formula, the solution never converges to a reasonable position. I don't know why.
Honestly speaking I'm not 100% sure if the solution above is reliable.
Comments
Post a Comment