This comment is written using Firefox.

Under Safari it is (at least for me) impossible to login.
Things start with an a priori satisfactory connection:





Now trying to send you a comment results in this pop-up window




And things keep going on this way.
Trying to asx a quastion from Safari results in an endless loop:





Clicking on Sign in reopen windows



The top right region of Mapleprimes indeed shows I'm not logged in after having  Sign in

@Math-dashti 

The file, with the output I asked you to produce, has too large a size.
No problem: just reduce the time range for instance replace

@Math-dashti 

PS: I'm mmcdara answering you from my professional account.

It's likely that the structure of dsol has changed betwwen Maple 2015 and Maple 2024.
To be sure of that can you simply insert this line just before the block  T := [seq(....)]

dsol;

and upload the worksheet with the lengthy ourput this command generates (do not execute the code beyond the block where T is assigned).

TIA

@Andiguys 

Optimisation package does not handle strict inequalities.
So it is impossible to impose, for instance, the condition w > Pu.
But writing w >= Pu+epsilon for a small value of epsilon should lead to the same solution, at least as long as constraints C1 or C2, or the objective function itself are not singular when epsilon -> 0⁺.

For instance both C1,  C2 and TP_2 are singular at upsilon=0 (so, I guess, your condition upsilon > 0).
So replacing upsilon > 0  by upsilon >= epsilon could lead to numerical problems.
Happily this does not seem to be the case here and using non strict inequalities with a "small"  epsilon value should indeed lead to the same solution than strict inequalities.

@salim-barzani 

You haave to produce images with the "quality" that the journal you target asks you to deliver

First select the journal you want to submit your paper to and read the corresponding guidelines.

For instance, guideline for the Journal of Fluid Mechanics can be found here with a special topic at parapraph 5.1 Figures.

All journals with an editing board and a review commitee editthose guidelines (for the others, for instance in rapidly evolving subjects like AI or Deep kearning, you can almost do what you want).

@JAMET 

If it so look at this Youtube video.
The first conic the speaker talks about is a parabola (your case).
Follow step by step what he does:

1. use LinearAlgebra:-Eigenvectors to diagonalize the Matrix(2$2, [9, -12, -12, 16]) associated to your quadratic form, let P the matrix of right eigenvectors, use the transformation <u, v> = P^(-1).<x, y> to align the parabola axis with u, (or v see below) ;
    
2. use Student:-Precalculus:-CompleteSquare to get the expression a*(u-p)^2 + b*v + c (or a'*(v-p')^2 + b'*u + c' depending on the order LinearAlgebra:-Eigenvectors returns the eigenvalues);
    
3. rewrite b*v + c into b*(v + c/b)  (or b'*(u + c'/b')) to get the translation aling u (v);
    
4. set U = u-p and V=v-c/b (or V = v-p' and U=u-c'/b' ) to get the reduced form a*U^2+b*V=0 (or a'*V^2+b'*U=0)
    

@dharr 

If I have some time I will try to design a more subtle algorithm.
 

@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 (X_i, Y_i, F(X_i, Y_i)). There is no F(X_i, Y_i) 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): (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)]) ); (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 ?

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