sand15

792 Reputation

11 Badges

9 years, 231 days

MaplePrimes Activity


These are replies submitted by sand15

@schachar 

You wrote " Please use the convolution of field (X, Y) by a 2D  gaussian for the data", but you did not provide any field.
You provided two series X and Y index by n (n=1..781).
A discrete 2D field F(X, Y) is a set of triples (Xi, Yi, F(Xi, Yi
)). There is no F(Xi, Yi) here.

@vv 

The epsilon was intended to be a workaround to the conditions r[i] > 0.
My intention was to write instead r[i] >= epsilon, that I changed to r[i]+epsilon >= 0 by carelessness.

Thanks for pointing out this stupid mistake.

@minhthien2016 

I got 1.905667204, so almost the same thing.

Do you know what algorithm your friend used (the (apparent?) symmetry with respect to the second diagonal is remarkable)?
Are solutions obtained by Asymptote for other values of n still symmetric with respect to a diagonal?
Those I got (see  mmcdara's reply) are symmetric for n from 1 to 9 only.

Meanwhile I edited my reply because the plots were not correct. You will se that the plots are now identical to the one you get with Asymptote.
 

By the way sand15 / mmcdara is one and the same person (provessional/personal login)

 

@acer 

I didn't check.
My intuition is that there should be equalities.

The situation would be different if you search the best way to pack N disks or same radius R into the unit square. In this case they may exist somewhere a free volume of size > Pi*r^2 where a last disk can be placed in an infinity ways, not even touching the other ones nor the boundary of the box.
But if you allow this last disk to have a radius R' > R to "fill the void", then there exist only one posibility to place it when R' reaches some maximum value.

 

@vv @Carl Love

EDITED
 

Some digging on the n=14 case.
I ran 3 sequences of 1000 optimizations all starting from a different initialpoint randomly chosen (unless if I am extremely unlucky all those initial points respect the constraints).

For these three sequences I kept the packing which gives the highest sum of radii (S).

I got three times the same value S=1.90567, the same sets of radii, and the packings are identical up to rotations by an angle multiple of Pi/2.

I belive these results partly answer @Carl Love's question "Can we confirm that the maxima are global?".
Indeed, I would infer that S=1.90567 is a global maximum.

As the two histograms provided seem to demonstrate, there could exist a large number of local maxima where NLPSolve is stucked.

 

restart

with(plots):
with(plottools):
with(Statistics):


DISKS WITHIN THE UNIT SQUARE

The packing shiuld be such that the sum of the radii is maximum

epsilon := 1e-8:

n := 14:
J := add(r[i], i=1..n):
P := [seq(<x[i], y[i]>, i=1..n)]:
cons_1 := seq(seq(add((P[i]-P[j])^~2) >= (r[i]+r[j])^2, j=i+1..n), i=1..n-1):
cons_2 := seq(op([P[i][1]-r[i] >= 0, P[i][1]+r[i] <= 1, P[i][2]-r[i] >= 0, P[i][2]+r[i] <= 1]), i=1..n):
cons_3 := seq(r[i] >= epsilon, i=1..n):


SoR : Sum of the Radii

Max_trials : maximum number of trials
                     Note: I used a while loop instead of a for loop because NLPSolve sometimes returns an error

Packing := proc(n, Max_trials)
  local SoR, SoR_max, error_found, error_count, rep, StartAt, opt, sor, error_type, Best_packing:

  SoR     := Vector(Max_trials):
  SoR_max := 0:

  error_found := false:
  error_count := 0:

  rep := 0:
  while rep < Max_trials do
    StartAt := {
                 seq(x[i]=rand(0. .. 1.)(), i=1..n),
                 seq(y[i]=rand(0. .. 1.)(), i=1..n),
                 seq(r[i]=rand(0. .. 1e-6)(), i=1..n)
               };
  
    try
      opt := Optimization:-NLPSolve(J, {cons_1, cons_2, cons_3}, maximize, iterationlimit=10^4, initialpoint=StartAt):
    catch:
      error_found := true
    end try:
  
    if `not`(error_found) then
      rep := rep+1:
      sor := evalf[6](eval(add(r[i], i=1..n), opt[2]));
      SoR[rep] := sor:
  
      DocumentTools:-Tabulate(
        [sprintf("%d", rep), sprintf("%1.6f", SoR_max)],
        'alignment'='left',
        'widthmode'='percentage',
        'width'=20
      ):
      
      if sor > SoR_max then
        SoR_max      := sor;
        Best_packing := opt[2];
      end if:
    else
      error_found := false:
      error_count := error_count+1:
      error_type  := lastexception:
    end if:
  end do:

  return  table([
            Errors = `if`(error_count = 0, None, [Type = error_type, Occurrences = error_count]),
            RadSum = SoR,
            Best   = Best_packing
          ])
end proc:

res := Packing(14, 1000):

res[Errors];

SoR := res[RadSum]:
Best_packing := res[Best]:

Best_packing_1 := copy(Best_packing):

[Type = (Optimization:-NLPSolve, "no improved point could be found"), Occurrences = 134]

(1)

FivePointSummary([SoR]):
Histogram([SoR], minbins=100)

 

interface(rtablesize=n+2):

`<,>`(
  `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),
  convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),
  `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)])
);

Matrix([[i, x[i], y[i], r[i]], [1, .464173791911759, .133714453381381, .133714453381381], [2, .833892835887959, .833892835887959, .166107164112041], [3, .650909465823001, .635654995533526, .103672662911138], [4, .166107164112041, .166107164112041, .166107164112041], [5, .866285546618619, .535826208088240, .133714453381381], [6, .715313251648918, .117921112194915, .117921112194915], [7, .915250635290483, 0.847493647095166e-1, 0.847493647095167e-1], [8, .519468470480418, .851206475514818, .148793524485182], [9, .882078887805085, .284686748351082, .117921112194915], [10, .623290519363489, .376709480636510, .156741595317294], [11, .364345004466473, .349090534176998, .103672662911138], [12, .148793524485182, .480531529519582, .148793524485182], [13, .415327133247299, .584672866752700, .137363045805254], [14, .186395365293329, .813604634706671, .186395365293329], [sum, ` `, ` `, 1.90566720529471]])

(2)

Histogram(
  (evalf[6]@rhs)~(select(has, Best_packing, r))
  , discrete=true
  , size=[1200, 100]
  , color=red
  , thickness=3
  , title="Radii"
);

Packing_1 := display(
  seq(
    disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)
    , i=1..n
  ),
  rectangle([0, 0], [1, 1], color="LightGray")
  , title=typeset(Sum(r[i], i=1..n) = evalf[6](eval(add(r[i], i=1..n), Best_packing)))
  , size=[500, 500]
  , scaling=constrained
):

 


RUN 2

res := Packing(14, 1000):

res[Errors]:

SoR := res[RadSum]:
Best_packing := res[Best]:

Best_packing_2 := copy(Best_packing):

FivePointSummary([SoR]):
Histogram([SoR], minbins=100):

`<,>`(
  `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),
  convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),
  `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)])
):

Histogram(
  (evalf[6]@rhs)~(select(has, Best_packing, r))
  , discrete=true
  , size=[1200, 100]
  , color=red
  , thickness=3
  , title="Radii"
);

Packing_2 := display(
  seq(
    disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)
    , i=1..n
  ),
  rectangle([0, 0], [1, 1], color="LightGray")
  , title=typeset(Sum(r[i], i=1..n) = eval(add(r[i], i=1..n), Best_packing))
  , size=[500, 500]
  , scaling=constrained
):

 


RUN 3

res := Packing(14, 1000):

res[Errors]:

SoR := res[RadSum]:
Best_packing := res[Best]:

Best_packing_3 := copy(Best_packing):

FivePointSummary([SoR]):
Histogram([SoR], minbins=100):

`<,>`(
  `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),
  convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),
  `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)])
):

Histogram(
  (evalf[6]@rhs)~(select(has, Best_packing, r))
  , discrete=true
  , size=[1200, 100]
  , color=red
  , thickness=3
  , title="Radii"
);

Packing_3 := display(
  seq(
    disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)
    , i=1..n
  ),
  rectangle([0, 0], [1, 1], color="LightGray")
  , title=typeset(Sum(r[i], i=1..n) = eval(add(r[i], i=1..n), Best_packing))
  , size=[500, 500]
  , scaling=constrained
):

 


The three best solutions found are identical after ad hoc rotations

DocumentTools:-Tabulate([
  Packing_1,
  translate(rotate(translate(Packing_2, -1/2, -1/2), -Pi/2), 1/2, 1/2),
  translate(rotate(translate(Packing_3, -1/2, -1/2), -Pi), 1/2, 1/2)
])

 

 

Download Disks_within_unit_square_14_edited.mw
 

 

@mehran rajabi 

Is the expression of f given once and for all or is it just a notional example that you provide?

Do you want Something_like_this.mw ?

restart

f := (u, v) -> u+v

proc (u, v) options operator, arrow; u+v end proc

(1)

rel0 := diff(y(x), x) = f(x, y(x)):
rel  := copy(rel0):

for n from 1 to 3 do
  rel  := rel, diff(lhs(rel0), x$n) = eval(expand(diff(rhs(rel0), x$n)), {rel});
end do:

print~([rel]):

   

diff(y(x), x) = x+y(x)

 

diff(diff(y(x), x), x) = 1+x+y(x)

 

diff(diff(diff(y(x), x), x), x) = 1+x+y(x)

 

diff(diff(diff(diff(y(x), x), x), x), x) = 1+x+y(x)

(2)
 

 

Here are A_few_generalizations.mw

Run Maple, open a new worksheet or document and execute this command

help(dsolve)

Read the help page and learn how to solve a differential equation with Maple.
Hint:

# y' = y+x means y'(x) = x+y(x) and is written diff(y(x), x) = x+y(x)
#
# At this point there is no need to use complex stuff like declaring that y is an alias for y(x) or
# using special tricks to us y' instead of diff(y(x), x)

dsolve(diff(y(x), x) = x+y(x))
                   y(x) = -x - 1 + exp(x) _C1

upload your worsheet by using the green arrow in the menu bar


UPDATED

This file contains the solutions of 3 well-posed problems:

  1. {f(0) = 0, (D(f))(0) = 0, ((D@@2)(f))(0) = a, theta(0) = 1, (D(theta))(0) = b}
  2. {f(0)=0, D(f)(0)=0, f(10)=1, theta(0)=1, theta(10)=0}
    
  3. {f(0)=0, D(f)(0)=0, D(f)(10)=1, theta(0) = 1, theta(10) = 0}
    

 

Make your choice.
And if none of these three problems corerspond to what you want, try to take inspiration of what is done in the attached file to solve your problem of interest.

Several_well-posed-problems_sand15.mw


@acer 

Do you think it's decent to improve your reputation based on work others have done?
You don't need to. You've already earned it.
Personally, I would have mentioned Kitonum's work in a simple reply, not an answer.

@janhardo 

As ChapGPT is meant to provide clear synthesis on a subject, the least one can hope for is that its answers are error free and understandable.
For instance

  • Can you explain me how ChatGPT passes from the 4th equality to the 5th one (in point 4)?
  • Can you explain me how ChatGPT finds that (point 7) pi3/3! = pi2/6?
  • More generally, ChatGPT does not shine for its clarity (Point 6 is quite obscure).


Better sources than CheatGPT exist  to understand how Euler solved the Basel problem.
Wiki (see @vv) is already better but nothing is worth the original  Euler (english translation)
 

@vv 

Thank you for this explanation (even if I only remember the analytical continuation in name, it's all too far away for me).
If I have a bit of time, I'll try to refresh my memory.

@Jesús Guillera 

By the way, did you browse the references I provided you or did you made your own research about convergence conditions of your serie?

To end my contribution: are you sure what you said in your last reply is correct?
 

kernelopts(version);
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895

F(-2);
                   infinity                 
                    -----                   
                     \              2     n 
                      )     po(1, n)  (-2)  
                     /     -----------------
                    -----          2        
                    n = 1     /1   \         
                           po| -, n |  (n + 1)
                              \2   /         

A truncated serie:

FT := (z, N) -> add(po(1,n)^2/po(1/2,n)^2*z^n/(n+1),n=1..N):
evalf(F(1/2)), evalf(FT(1/2, 100)); 
                    2.254576561, 2.254576561

 evalf(F(-2)) , evalf(FT(-2, 10)), evalf(FT(-2, 50)), evalf(FT(-2, 250)); 
                                          15                75
 -1.076461643, 2002.463620, 2.323658088 10  , 3.777981621 10  

Convergence doesn't seem to hold at z=-2...

Look at this: Download Convergence_criteria.mw

@Jesús Guillera 

I'm not a specialist of hypergeometric series, but if you're right then this is a Maple bug.

F(2)
                            infinity


But on the other side, if you "think that F(2) is convergent" why do you declare assume(abs(z)<1) ?


This communication (unfortunately in French) should put an end to this dicussion, mid of page 5 "
Le rayon de convergence est évidemment l'unité"  which translate as "The convergence radius is obviously the unity" meaning the serie is convergent if abs(z) < 1 (abs(z) <= 1 ?).

An hypergeometric serie

 F := (z, a, b) ->sum(po(a,n)^2/po(b,n)^2*z^n/(n+1),n=1..infinity)

is convergent for any real z if at least one of a or b is a strictly negative integer (If I'm not mistaken, because the associated Pochammer symbol cancels out beyond a certain value of n).
But you are not in this situation with  a=1 and b=1/2 ...

See also Rivoal (still in french) from which this excerpt is taken:



The highlight text translate as "with sometimes a convergence on the unit circle"


So (once again I'm not a specialist), I advice you to verify if your hypergeometric serie converges for z > 1.

@C_R 

I vote up.

In fact I spent a little time on this question where I followed a slightly different path.
But I stpped digging after some time

1 2 3 4 5 6 7 Last Page 1 of 25