MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

Mathematica program to obtain a bounding function for a set of data points.

  • To: mathgroup at smc.vnet.net
  • Subject: [mg43248] Mathematica program to obtain a bounding function for a set of data points.
  • From: gilmar.rodriguez at nwfwmd.state.fl.us (Gilmar Rodríguez Pierluissi)
  • Date: Sat, 23 Aug 2003 08:08:23 -0400 (EDT)
  • Sender: owner-wri-mathgroup at wolfram.com

The following is to help clarify my question about how to build a
Mathematica
program to obtain a bounding function for a set of data points (please
refer to:
http://forums.wolfram.com/mathgroup/archive/2003/Jul/msg00496.html).
I haven't received a response so far, and I figure that the reason for
this is that I need to elaborate my question better.

(** Preliminary (but, necessary)information :**)

Remark: Q_i  below means "Q subscript i", etc.

Let N be an integer greater or equal to 4. Assume the existance of
primes P_i
in the interval [2, N/2] in X, and primes Q_i  = (N - P_i) in the
interval
[N/2, N - 2] in Y, with P_i <= Q_i such that N = P_i + Q_i. Then the
set of points(P_i,Q_i) is called a Prime Partition of N.

The Minimal Goldbach Decomposition Point ("MGDP" for short) 
is a point (P,Q) among the points (P_i,Q_i) such that the
perpendicular distance
between the points (P_i,Q_i) (located on the line Y = -X + N) and the
point (N/2,N/2)(located on the line Y=X) is the shortest distance,
i.e.

perpdist[(P,Q), (N/2,N/2)]=Min[Sqrt[(P_i - N/2)^2+(Q_i- N/2)^2]]

The following web pages help to clarify the concepts considered so far
:

http :    www.gilmarlily.net goldbach gcstated.htm

http :    www.gilmarlily.net goldbach geoview.htm

http :    www.gilmarlily.net goldbach mgdplot1.htm

(** End of preliminary stuff. **)


First step; the following module gives the Minimal Goldbach
Decomposition for an integer n >= 4.

In[1]:

MGDP[n_] := Module[{p, q},
    {m = n 2; 
      If[(m \[Element] Primes), {(p = m), (q = m)}, {k = PrimePi[m]; 
          Do[If[(n - Prime[i]) \[Element] Primes, {hit = i, Break[]}],
{i, k,
              1, -1}], p = Prime[hit], q = (n - p)}]}; {p, q}]

For example; 

In[2]:

MGDP[100]

Out[2]:
{47, 53}

Second step; we build a Module to calculate the perpendicular distance
between
the MGDP and the point (n/2, n/2). For each n, the Module returns a
value
{n, perpendicular distance}.

In[3]:

perpdist[n_] := Module[{p, q, d},
    {m = n 2; dist[{x1_, y1_}, {x2_, y2_}] := Sqrt[(x2 - x1)^2 + (y2 -
y1)^2];
       If[(m \[Element] Primes), {(p = m), (q = m), (d = 0)}, {k =
PrimePi[m];
           Do[If[(n - Prime[i]) \[Element] Primes, {advpt = i,
Break[]}], {i,
              k, 1, -1}], p = Prime[advpt], q = (n - p), 
          d = dist[{p, q}, {m, m}]}]}; {n, d}]

For example;

In[4]:

perpdist[100]

Out[4]:

{100, 3Sqrt[2]}

Third step; for each n between 4 and 1 Million, we calculate the
perpendicular
distance between the MGDP of n, and the point (n/2,n/2).  (The
calculation of the set A might take a little time):

In[5]:

A = Table[perpdist[n], {n, 4, 10^6, 2}];

Fourth step; we seek to find a Bounding Function for the set A above.

To start, we first build a module called SMaxima[]. This module
compares
each pair {(a_x,a_y), (a_(x+1),a_(y+1))} in A, and appends
(a_(x+1),a_(y+1)) to a set M, if a_y < a_(y+1).

In[6]:

SMaxima[B_] := 
  Module[{M}, {M = {}; AppendTo[M, B[[1]]]; {Xg, Yg} = B[[1]]; 
      Do[If[Yg < B[[i + 1]][[2]], {{Xg, Yg} = B[[i + 1]], 
            AppendTo[M, B[[i + 1]]]}], {i, 1, Length[B] - 1}]}; M]

Applying SMaxima[] to A gives us :

In[7]:

M = SMaxima[A]

Out[7]:

{{4, 0}, {8,  Sqrt[2]}, {16, 3Sqrt[2]}, {44, 9Sqrt[2]}, {92,
15Sqrt[2]},
 {242, 18Sqrt[2]}, {256, 21Sqrt[2]}, {272, 27Sqrt[2]}, {292,
33Sqrt[2]},
 {476, 39Sqrt[2]}, {530, 42Sqrt[2]}, {572, 45Sqrt[2]}, {682,
48Sqrt[2]},
 {688, 75Sqrt[2]}, {1052, 87Sqrt[2]}, {1808, 93Sqrt[2]}, {2228,
117Sqrt[2]},
 {3382, 120Sqrt[2]}, {3472, 135Sqrt[2]}, {3502, 138Sqrt[2]},
 {3562, 168Sqrt[2]}, {4952, 183Sqrt[2]}, {6194, 210Sqrt[2]},
 {7102, 228Sqrt[2]}, {10262, 300Sqrt[2]}, {17008, 333Sqrt[2]},
 {20684, 369Sqrt[2]}, {37052, 393Sqrt[2]}, {45128, 453Sqrt[2]},
 {49552, 525Sqrt[2]}, {80144, 621Sqrt[2]}, {137414, 720Sqrt[2]},
 {251806, 810Sqrt[2]}, {349826, 846Sqrt[2]}, {362534, 1086Sqrt[2]},
 {742856, 1281Sqrt[2]}}

Plot M using the command:

In[8]:

plt1 = ListPlot[M, PlotJoined -> True, PlotRange -> All, PlotStyle ->
Hue[.6]]

Out[8]: 

You get a plot.

Fifth step; the set M becomes instrumental in finding a Bounding
Function for A :

In[9]:

<< Statistics`NonlinearFit`

We want our bounding function to pass through the point in M with the
highest
y - coordinate, so we do :

In[10]:

Off[General::"spell1"]

In[11]:

data = {{M[[Length[M]]][[1]], M[[Length[M]]][[2]]}}

Out[11]:

{{742856, 1281Sqrt[2]}}

In[12]:

BF[N_] = NonlinearFit[data, Alpha/(Exp[2/ Log[N]] - 1), {N}, {Alpha}]

Out[12]:

(288.8656741801377)/(-1 + Exp[2/Log[N])

We now build a set called "Bounding Set" (i.e. "BoundSet" for short)
where the points have the same x - coordinates as the points in M, 
but the y - coordinates are calculated using our Bounding Function
above :

In[13]:

BoundSet = Table[{M[[i]][[1]], BF[M[[i]][[1]]]}, {i, 1, Length[M]}]

Out[13]:

{{4, 89.3744}, {8, 178.71}, {16, 273.236}, {44, 414.792}, {92,
519.275},
 {242,657.102}, {256, 665.136}, {272, 673.799}, {292, 683.94},
 {476, 753.851}, {530, 769.239}, {572, 780.162}, {682, 805.363},
 {688, 806.618}, {1052, 867.505}, {1808, 945.222}, {2228, 975.218},
 {3382, 1035.18}, {3472, 1038.95}, {3502, 1040.19}, {3562, 1042.63},
 {4952, 1089.99}, {6194, 1122.17}, {7102, 1141.84}, {10262, 1194.79},
 {17008, 1267.49}, {20684, 1295.65}, {37052, 1379.59}, {45128,
1407.98},
 {49552, 1421.45}, {80144, 1490.7}, {137414, 1568.38}, {251806,
1655.66},
 {349826, 1703.05}, {362534, 1708.19}, {742856, 1811.61}}

We now plot the points in BoundSet :

In[14]:

plt2 = ListPlot[BoundSet, PlotJoined -> True, PlotRange -> All, 
    PlotStyle -> Hue[.4]]

Out[14]:

You get a second plot.

Finally, we plot the sets M and BoundSet together :

In[15]:

Show[plt1, plt2]

Out[15]:

You get a third plot.

The point that I'm really trying to make via this example is : 
it would be nice if Wolfram Research could combine the "SMaxima[]
technique" above, together with the NonlinearFit program to come up
with a
new program that calculates the Bounding Function for an arbitrary 
set of
points similar to the set A above.  Thank you for your attention and
your help!


  • Prev by Date: Working with Arbitrary Precision and NDSOLVE
  • Next by Date: Re: Re: New version, new bugs
  • Previous by thread: Working with Arbitrary Precision and NDSOLVE
  • Next by thread: Mesh in one direction only?