Unable to Solve Area Using Paramteric Equations derived from The Numeric Soution of a Differential Equation

I'm creating a test that measures the distribution of closed shape $\mathcal{C}$. This requires taking the average ray length (from $(u,v)$ in $\mathcal{C}$ to the boundary) using ray length per radian measure.

Suppose $\mathcal{C}$ can be described as $(a_1(t),a_2(t))$. If $\overline{r}(u,v)$ is the average radius, the distrbution number is

$$\phi(u,v)=\frac{1}{\overline{r}(u,v)}\int\limits_{t_1}^{t_2}\left( \left|\sqrt{\left(a_1(t)-u\right)^2+\left(a_2(t)-v\right)^2}-\overline{r}(u,v)\right|\frac{d\theta}{dt}\left(\frac{1}{2\pi}\tan^{-1}\left(\frac{a_2(t)-u}{a_1(t)-v}\right)\right) \right)dt$$

Where $t_1<t<t_2$ parametrizes the entire region

There is no clear way of determining where to measure $(u,v)$. Hence I took the average of $\phi(u,v)$ for all $(u,v)$ in $\mathcal{C}$. This leads to a difficult formula.

$$\omega=\frac{1}{\text{Area}\left(\mathcal{C}\right)} \iint\limits_{C} \phi(u,v) \ du \ dv$$

In this answer, I can parametrize any implicitly defined $\mathcal{C}$ using arclength. In this example, $\mathcal{C}$ is $\left(x^2 + y^2 +\sin(4 x) + \sin(4 y) - 4\right)^2 + .415 x=.4$

expr = (x^2 + y^2 + Sin[4 x] + Sin[4 y] - 4)^2 + .415 x - .4;
exprt = expr /. {x -> x[t], y -> y[t]};
odes = {D[exprt, t], x'[t]^2 + y'[t]^2 - 1}
starty = y /. FindRoot[(expr /. x -> 0) == 0, {y, 2/9}]
Quiet[soln =
NDSolveValue[
Join[Thread[odes == 0], {x == 0, y == starty}], {x[t],
y[t]}, {t, 0, 25}];]
{solnX, solnY} = soln;
ParametricPlot[soln, {t, 0, 25}]
Plot[(solnX - 0)^2 + (solnY - starty)^2, {t, 0, 25}]
{min, val} = NMinimize[(solnX - 0)^2 + (solnY - starty)^2, {t, 20, 25}]
stopt = t /. val;
r = Compile[{{x, _Real}, {y, _Real}},
NIntegrate[
Sqrt[(solnX - x)^2 + (solnY - y)^2]*
D[ArcTan[solnX - x, solnY - y], t], {t, 0, stopt},
MaxRecursion -> 100, Method -> "Trapezoidal"]/(2*Pi)]

s = Compile[{{x, _Real}, {y, _Real}},
NIntegrate[
Abs[Sqrt[(solnX - x)^2 + (solnY - y)^2] - r[x, y]]*
D[ArcTan[solnX - x, solnY - y], t], {t, 0, stopt},
MaxRecursion -> 100, Method -> "Trapezoidal"]/(2*Pi)]

The distrbution number is

Phi = Compile[{{x, _Real}, {y, _Real}}, s[x, y]/r[x, y]]

The problems is I'm unable to use the code in the following ways.

1. I can't extract values from solnX and solnY. I tried solnX, solnX[] and sol[{{0.,25.}}]. (Note these are not "typical" interpolation formulas, they're composed of two parts: the domain output and t.) This was fixed in the comments below.

1. I cant graph Graphing Phi[x,y] takes too long.

2. I can't solve for the area using Integrate[D[solnX,t]*solY,{t,0,stopt}]

Without these, I cannot compute $\omega$.

• have you tried to make this work without monkeying around with Compile ? – george2079 Jan 26 '18 at 14:49
• @george2079 I tried without Compile, it didn't work. – Arbuja Jan 26 '18 at 15:05
• @george2079 The main problem is with SolnX and SolnY – Arbuja Jan 26 '18 at 15:14
• do like this: soln[t_] = NDSolveValue... , then soln -> {-0.54648, 2.12265} or you could do {solnX[t_], solnY[t_]} = soln with soln the way you have it. – george2079 Jan 26 '18 at 15:17
• Thanks, it worked. – Arbuja Jan 26 '18 at 15:26