mmcdara

7891 Reputation

22 Badges

9 years, 56 days

MaplePrimes Activity


These are replies submitted by mmcdara

@sursumCorda 

Thanks for this information.
My version of Maple is too old to run correctly this code. I'll look at it when back at the office.

Here is a procedure aimed to find the tangent points between the graph of a curve of expression f(x,y)=0 and circles of equation 
x^2+y^2-R^2=0.
This procedure returns a list of three elements:

  1. the triples (x, y, R)
  2. the minimum value of R (tangency points closest to the origin)
  3. the coordinates of the tangency points closest to the origin

The procedure needs to be improved for it gives wrong result with function 
f(x, y)  = x^3 + y^2 - x + 2/(3*sqrt(3))  
(an isolated point [x=1/sqrt(3), y=0] is detected which is obviously not on the graph of f(x, y))

restart;

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

with(plots):

ClosestPoint := proc(F)
  local f, g, Df, Dg, sys, sol, xR, Y, xyR, i, j:
  local Rmin, MPP:

  if F::`=` then
    f := (lhs-rhs)(F)
  else
    f := F:
  end if:
  g := x^2+y^2-R^2;

  Df := diff~(f, [x, y]);
  Df := Df /~ sqrt(Df[1]^2+Df[2]^2);

  g  := x^2+y^2-R^2;
  Dg := diff~(g, [x, y]);
  Dg := Dg /~ sqrt(Dg[1]^2+Dg[2]^2);

  sys := {(Df =~ Dg)[], f=g};
  sol := solve(sys, [x, y, R])[];

  # contact points with vertical tangent
  {solve({Df[2]=0}, [x, y])[]};
  allvalues(%);
  map(s -> solve(eval(f, s), x), %);
  remove(has, simplify~(%), I);
  sol := sol, map(s -> [x=s, R=abs(s), y=0], %)[];


  sol := map(allvalues, [sol]):

  xR  := map(s -> if (`not`@has)(eval(R, s), x) then remove(has, s, y) end if, sol);
  Y   := map(s -> [solve(eval(g, s), y)], xR);
  xyR := NULL:
  for i from 1 to numelems(xR) do
    for j from 1 to numelems(Y[i]) do
      xyR := xyR, [ xR[i][], y=Y[i][j] ]:
    end do:
  end do:
  xyR := remove(has, {xyR}, I);
  xyR := map(s -> if evalf(eval(R, s)) > 0 then s end if, xyR);

  Rmin := min(eval~(R, xyR));
  MPP  := select(has, xyR, (R=Rmin));

  return [ xyR, Rmin, MPP ]
end proc:

C := ClosestPoint(x^3+y^2-x=1/3)

[{[x = 1, R = (2/3)*3^(1/2), y = -(1/3)*3^(1/2)], [x = 1, R = (2/3)*3^(1/2), y = (1/3)*3^(1/2)], [x = -1/3, R = (2/9)*3^(1/2), y = -(1/9)*3^(1/2)], [x = -1/3, R = (2/9)*3^(1/2), y = (1/9)*3^(1/2)], [x = -sin((1/18)*Pi)-(1/3)*3^(1/2)*cos((1/18)*Pi), R = sin((1/18)*Pi)+(1/3)*3^(1/2)*cos((1/18)*Pi), y = 0], [x = sin((1/18)*Pi)-(1/3)*3^(1/2)*cos((1/18)*Pi), R = -sin((1/18)*Pi)+(1/3)*3^(1/2)*cos((1/18)*Pi), y = 0]}, (2/9)*3^(1/2), {[x = -1/3, R = (2/9)*3^(1/2), y = -(1/9)*3^(1/2)], [x = -1/3, R = (2/9)*3^(1/2), y = (1/9)*3^(1/2)]}]

(2)

display(
  implicitplot(x^3+y^2-x = 1/3, x=-2..2, y=-2..2, grid=[100, 100], color=red)
  , seq(plottools:-circle([0, 0], eval(R, c), color=blue), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=box, symbolsize=20, color=magenta), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=diagonalcross, symbolsize=20, color=magenta), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=box, symbolsize=20, color=green), c in C[3])
  , seq(plot([eval([x, y], c)], style=point, symbol=diagonalcross, symbolsize=20, color=green), c in C[3])
  , scaling=constrained
  , view=[-1.1..1.5, -1..1]
  , title="Green points are those closest to the origin"
);

 

C := ClosestPoint(x^3+y^2-x = -2/(3*sqrt(3)))

[{[x = -(2/3)*3^(1/2), R = (2/3)*3^(1/2), y = 0], [x = (1/3)*3^(1/2), R = (1/3)*3^(1/2), y = 0]}, (1/3)*3^(1/2), {[x = (1/3)*3^(1/2), R = (1/3)*3^(1/2), y = 0]}]

(3)

display(
  implicitplot(x^3+y^2-x = -2/(3*sqrt(3)), x=-2..2, y=-2..2, grid=[100, 100])
  , seq(plottools:-circle([0, 0], eval(R, c), color=blue), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=box, symbolsize=20, color=magenta), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=diagonalcross, symbolsize=20, color=magenta), c in C[1])
  , seq(plot([eval([x, y], c)], style=point, symbol=box, symbolsize=20, color=green), c in C[3])
  , seq(plot([eval([x, y], c)], style=point, symbol=diagonalcross, symbolsize=20, color=green), c in C[3])
  , scaling=constrained
  , view=[-2..2, -2..2]
  , title="Green points are those closest to the origin"
);

 

# Why is the result wrong ?

  sol := NULL:

  f  := x^3+y^2-x+2/(3*sqrt(3));
  Df := diff~(f, [x, y]);

  # contact points with vertical tangent

  {solve({Df[2]=0}, [x, y])[]};
  allvalues(%);
  map(s -> solve(eval(f, s), x), %);
  remove(has, simplify~(%), I);
  sol := map(s -> [x=s, R=abs(s), y=0], %)[];

x^3+y^2-x+(2/9)*3^(1/2)

 

[3*x^2-1, 2*y]

 

{[x = x, y = 0]}

 

{[x = x, y = 0]}

 

{-(2/3)*3^(1/2), (1/3)*3^(1/2)}

 

{-(2/3)*3^(1/2), (1/3)*3^(1/2)}

 

[x = -(2/3)*3^(1/2), R = (2/3)*3^(1/2), y = 0], [x = (1/3)*3^(1/2), R = (1/3)*3^(1/2), y = 0]

(4)

# Point [x = (1/3)*sqrt(3), R = (1/3)*sqrt(3), y = 0] is the wrong

eval(f, [x = (1/3)*sqrt(3), R = (1/3)*sqrt(3), y = 0])

0

(5)

 

Download ClosestPoint.mw

As the second function f(x, y) has a null genus, it admits a parameterized representation [x(t), y(t)].
It appears that this approach here returns the correct answer:

f  := x^3+y^2-x + 2/(3*sqrt(3)); 
                      3    2       2  (1/2)
                     x  + y  - x + - 3     
                                   9       
Gf := algcurves:-genus(f, x, y); 
                               0
if Gf = 0 then   
  A   := algcurves:-parametrization(f,x,y,t);  
  C   := R *~ [(1-s^2)/(1+s^2), 2*t/(1+s^2)];   
  sys := [(diff(A, t) =~ diff(C, t))[], (A =~ C)[]]:   
  sol := solve(sys, [s, t, R]);   
  MPP := [x, y] =~ eval(A, sol[]); 
end if: 

MPP; 
      [       1  /   (1/2)    \ /          (1/2)\       ]
      [x = - --- \2 3      + 9/ \-12 + 18 3     /, y = 0]
      [      207                                        ]

expand(simplify(%)); 
                  [x = -(2/3)*sqrt(3), y = 0]
evalf(%);
                  [x = -1.154700539, y = 0.]

This observation could be used to improve the robusteness of procedure ClosestPoint.


just a simpler version of your code (I was completely unaware of Draghilev's method)
Draghilev_mmcdara.mw

@acer 

Thanks for your reply.
Couldn't it be possible that Mapleprimes offers a kind of "log file" to trace such movements?
We never know with certainty if the disappearance of a question comes from a doubloon with another one, a conversion in post, or a suppression by the author himself

@sursumCorda 

Hi, what happened to your last question about extrema (tryHarder.mws)?
Did you delete or move it elsewhere, or did someone else took the initiative?

It always f.....g sucks to spend time preparing an answer and then realize you would have been better off taking a leak. 
(My apologies to those who felt shocked by these words)

It works functionally well with Maple 2015.2... but  the result is wrong whatever the version

Download Closest_Point.mw
 

For your information, your problem is a classic of structural reliability analysis.
The "standards" algorithms used to solve this kind of problems are:

  • The Abdo-Rackwitz (AR) algorithm.
    T. ABDO, R. RACKWITZ, A new beta-point algorithm for large time-invariant and time-variant reliability problems, [in:] A. DER KIUREGHIAN and P. THOFT-CHRISTENSEN [Eds.], Reliability and Optimization of Structural Systems '90 Proc. 3rd WG 7.5 IFIP Conf. Berkeley 26–28 March 1990, 1–12, Berlin 1991.
    (Didn(t find an open source)
  • The Hassofer-Lind-Rackwitz-Fiessler (HLRF) algorithm (a variant of SQP)
    https://downloads.hindawi.com/journals/mpe/2017/4317670.pdf
    https://www.researchgate.net/publication/283997508_New_optimization_algorithms_for_structural_reliability

@dharr 

Thanks for the link

@dharr 

Taking up your original idea (I voted up), here is a little bit faster way to get the solution

with(inttrans):

Lsys := eval(laplace(sys, t, s), ics): 
F    := convert(remove(has, indets(sys, function), diff), list):
F2L  := map(u -> op(0, u)(s), F):
L    := laplace~(F, t, s): 
Asys := eval(Lsys, L =~ F2L): 
Asol := solve(Asys, F2L)[]: 
Fsol := F =~ invlaplace~(rhs~(Asol), s, t):

Download dsolve2_mmcdara.mw
 

@dharr 

Hi, 
A passing remark: have you noticed how combinat:-cartprod is slow compared to your simple double loop ?

@sursumCorda 

I don't think your comparison is very reliable.
To enhance the credibility of CodeTools:-Usage outputs it is advisable touse the option iterations.
For a single iteration the cpu times vary from 8ms to 12ms (note this is muchsmaller than the 47ms you got).
For 100 iterations the same times vary from 9.88ms to 10.44ms.
One can say the mean cpu time is a little bit below 10ms.

It also appears that i||j is faster than cat(i, j) (third test).

It's also interesting to notice that the time used to construct the cartesian product of the two lists is even substantially smaller (fourth test)

restart;

interface(version);
interface(displayprecision=4):

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

T:= NULL:
for i from 1 to 1000 do
T := T, CodeTools:-Usage([seq](seq(cat(i,j),j in [$("a".."z"),$("A".."Z"),$("A".."Z"),$("a".."z"),$("A".."Z")]),i in [$("A".."Z"),$("a".."z"),$("a".."z"),$("A".."Z"),$("a".."z")]), iterations=1, output=cputime, quiet):
end do:

Statistics:-Histogram([T], discrete=true, thickness=10)

 

T:= NULL:
for i from 1 to 1000 do
T := T, CodeTools:-Usage([seq](seq(cat(i,j),j in [$("a".."z"),$("A".."Z"),$("A".."Z"),$("a".."z"),$("A".."Z")]),i in [$("A".."Z"),$("a".."z"),$("a".."z"),$("A".."Z"),$("a".."z")]), iterations=10, output=cputime, quiet):
end do:


Statistics:-Histogram([T], discrete=true, thickness=8)

 

T:= NULL:
for i from 1 to 1000 do
T := T, CodeTools:-Usage([seq](seq(i||j,j in [$("a".."z"),$("A".."Z"),$("A".."Z"),$("a".."z"),$("A".."Z")]),i in [$("A".."Z"),$("a".."z"),$("a".."z"),$("A".."Z"),$("a".."z")]), iterations=10, output=cputime, quiet):
end do:

Statistics:-Histogram([T], discrete=true, thickness=8)

 

L1 := [$("a".."z"),$("A".."Z"),$("A".."Z"),$("a".."z"),$("A".."Z")]:
L2 := [$("A".."Z"),$("a".."z"),$("a".."z"),$("A".."Z"),$("a".."z")]:

T:= NULL:
for i from 1 to 1000 do
T := T, CodeTools:-Usage([seq](seq(i||j, j in L1),i in L2), iterations=10, output=cputime, quiet):
end do:

Statistics:-Histogram([T], discrete=true, thickness=8)

 

N1 := numelems(L1):
N2 := numelems(L2):
A1 := Matrix(N1, N2, (i, j) -> L1[i]):
A2 := Matrix(N1, N2, (i, j) -> L2[j]):


T:= NULL:
for i from 1 to 1000 do
T := T, CodeTools:-Usage(cat~(A1, A2), iterations=10, output=cputime, quiet):
end do:

Statistics:-Histogram([T], discrete=true, thickness=8)

 

 

Download compareTime_mmcdara.mw

Unfortunately I can't make your code work in 2D math to compare the cpu times (note that the memory used is another important criterion for comparisons) :

interface(version); s1 := CodeTools:-Usage([cat("A" .. "Z"), cat("a" .. "z"), cat("a" .. "z"), cat("A" .. "Z"), cat("a" .. "z") || '`$`("a" .. "z"), `$`("A" .. "Z"), `$`("A" .. "Z"), `$`("a" .. "z"), `$`("A" .. "Z")'])

Error, unexpected single forward quote

 

``

Download 2Dmath.mw

In any case, I do prefer @dharr's method

@sursumCorda 

You hit the point: SAGE is free and extremely powerful thanks to the large number products its distribution includes (https://doc.sagemath.org/html/en/reference/spkg/)

About 15 years ago the main general purpose CAS used in French Universities (math) was Maple.
 I remember when Wolfram's commercials came presenting the last Mathematica release everyone was very critic: the show was all around the visualization, interactive features and a lot of other things mathematicians considered has "frills for engineers not for serious people". And they have kept their confidence in Maple.

But Maplesoft politics has changed and, still in France, when a Maplesoft's team comes to present Maple to mathematicians, they don't talk much about math but mostly about everything else, just as Wolfram's team did 10 years ago.
So there has been a slow but steady shift towards something less eye-catching, more focused on the things that these "serious people" consider essential... and less expensive
Today a lot of math's labs have turned towards SAGE, and some friends of mine, yet major Maple defendants have crossed the rubicon.

I suspect it's must be hard for Maple to find the right place between Mathematica and SAGE.

@sursumCorda 

Yes, R can do that.

Your initial question is interesting in that it points to a general problem: can Maple (or any other produc) be the best in any field?
And the answer is obviously not.

In your casse Mathematica and Matlab both produce better drawings of the network than Maple does.
But it's clear that Maple does other things far better than Matlab (which is not the same kind of product) and even Mathematica.
Trying to obtain with Maple what another product gives you without pain can be a huge task. As I use to use Maple and R, I have already recoded a lot of R-algorithms in Maple language, but it's worth it only if the price to pay for the development is not too big.

I'm more in favor of an hybrid strategy: if R (or Matlab, Python, ...) is "better" than Maple for doing some task, then use this product, instead use Maple.
As an example I keep using the Statistics  Maple's package to do formal computations or "simple" statistics, but switch on R do do complex statistical studies. But, as R itself has only limited capabilities on symbolic calculus (just as Matlab), I often swith from R to Maple.
Just trying to use the best tool for the task to do.

Information passing through files is not a problem. Maple can even delegate these tasks to R (or a lot of other products) through the ssystem/system mechanism.

So the question "I've obtained this result with ..., can Maple gets it" is a complex one.
I believe it is sometimes a legitimate question when the fiels of the products overlap (Mathematica solves this equation, can Maple do it?) but a pointless question if they dont (I solved this CFD problem with Fluent, can Maple do it?)
The pretty drawing of graphs or networks is not the strong point of Maple... but its GraphTheory package remains excellent.

@zenterix 

Here is how your file is displayed in 1D mode:

w := VectorCalculus:-RootedVector(root = [-2, 1], `<,>`(1, 2))

Vector[column](%id = 36893488152235632508)

(1)

VectorCalculus:-PlotVector(w, scaling = constrained)

 

VectorCalculus:-PlotVector(w)

 

NULL

Download PlotRootedVector.mw

Correct plots are displayed when this file is executed, but I can't find any plot command.

@zenterix 

I never use 2D math.
I don't consider it as reliable enough for me to use it (I did, years ago, and faced a lot of big problems, so I forgot about it).

But I kindly decided to make an exception :-)

with(VectorCalculus):

w := RootedVector(root = [-2, 1], `<,>`(1, 2)):

PlotVector(w, scaling = constrained)

 

``

Download 2D_math.mw

Then it seems to work also in 2D math

This is an example of what you can get with R (more customization is possible)
network.html

The document is interactive: use the green buttons or move the nodes

First 59 60 61 62 63 64 65 Last Page 61 of 154