@dharr I adjusted your Answer slightly, and inlined images (right-click exporting the plots) so that it would be certain that Mapleprime's own rendering wasn't interfering.  I hope that's OK.

@nm You have G as  x*sin(Pi*x)  instead of  x*sin(Pi*y) .

But even if your second result were to attain, it has x restricted to integers and hence not free.

@dondomingo But you won't know home many solutions there are and how many names to use, in general. I don't think that is a good methodology.

I agree, nobody wants to have to re-type or copy&paste from output. It's not good programming to have to do so.

Just assign the full list of solutions to a single name, let's call it Sol, say.

Then you can get at the first solution by indexing, as Sol[1]. And then you evaluate any formula involving x and y using 2-argument eval (as showed in my Answer above). You could even access just the simplest formula x, or y.

Eg,
   eval( x, Sol[1] );
   eval( y, Sol[1] );
   nops(Sol);    # how many solutions there are
   eval(  x^2+sin(y-x),  Sol[3] );
etc.

@alijp I don't understand why you didn't show that when you initially asked the Question.

Now you've introduced the `i` parameter when you already had two (r and phi). So you don't need `i`, now.

Is your problem with the differentiation at symbolic order n?

(I have marked your Question as being for Maple 17.)

I could point out that a function might by undefined on an interval between discontinuities (including a pair returned by discont). And iscont could return true for the test on that open interval. This might matter to you.

@Anthrazit Your string collection omitted the opening and closing square brackets, needed so that after parsing a list would be passed to Join.

joinStr_ac.mw

This is an example of what I meant about confusing multiple concatenation, since all the little fiddly strings are wrapped in their own quotes. I think that an sprintf format is more legible and less error-prone.

@Thomas Richard Thanks for the detail.

I mentioned those two particular versions because I believe those are the two that the OP uses.

@Carl Love For the sake of completeness, I'll mention that the collapsed calling sequence for a multiple integral can be more than just a syntactic convenience (versus the iterated calling sequence).

When computing the value, the ranges of the outer variables can be used during the computation of the inner integrals, via assumptions on the relevant names of the outer variables of integration. It turns out that this can quite often be useful, and the collapsed form is used in some VectorCalculus and MultivariateCalculus in Maple 2020.

@mmcdara I think that I didn't understand your original goal. Based on the original code you gave (and with those R commands slightly in mind) I had the impression that you were trying to get a kernel-density-plot based on given Samples, rather than from symbolic RVs.

Hence the suggestions to try surfdata with dimension=2, and direct writing to an Image format, etc. I didn't mention plots:-listdensityplot, which has a possibly useful smooth=false choice but lacks the convenience of surfdata to specify the ranges for the axes. I was working on the presumption that one would be starting from (only) samples of unknown distributions. I suppose there may be room for both goals.

It might be interesting to tune up your scheme for using PDF and the symbolid RVs. As is, each call to KD2D calls the K procedure N times, and each K call generates its own procedure as return value, and each call to those returns call the f procedure twice, and each call to the f procedure generates its own procedure as return value. Alas I don't have any time for that fun these days.

I notice that you bumped down your Sample size to N=10^3, perhaps for time cost and perhaps also because at N=10^4 the points occlude the filled color gradient of the density plot. And they collide so much that the visual gradient is lost. (BTW, you could switch the order of plots in the display, so that at least the black lines of the contours lay in front of the possibly dense scatter-point rendering.)

So here is a quick little thing. In the case of a large Sample there is some visual benefit to plotting a scattering using points smaller than symbolsize=1. By replacing the points with very thin and short lines we can get a more visually distinguishable mist.

>  restart: >  with(plots): with(Statistics): >  f := (m, s) -> unapply(PDF(Normal(m, s), t), t): K := (mx, my, sx, sy) -> (x, y) -> f(mx, sx)(x) * f(my, sy)(y): >  N := 10^3: A := Sample(RandomVariable(BetaDistribution(2, 3)), N*60): B := Sample(RandomVariable(BetaDistribution(3, 2)), N*60): h := 0.07: KD2D := (x, y) -> add(K(A[n], B[n], h, h)(x, y), n=1..N)/N: >  CP := contourplot(KD2D(x, y), x=0..1, y=0..1, view=[0..1, 0..1],                   filled=true, coloring=["MediumBlue", "White"],                   contourlabels=false, thickness=0): >  SP := ScatterPlot(A, B, symbol=point, color="Gold", symbolsize=1): >  dat := op([1,1],SP): newSP := subsop([1,1]=Matrix(3*op([1,1],dat),op([1,2],dat),                              (i,j)->`if`(irem(i,3)=1,dat[(i+2)/3,j],                                          `if`(irem(i,3)=2,1e-3+dat[(i+1)/3,j],                                               Float(undefined)))),                 subsindets(SP,specfunc(anything,STYLE),                            ()->[STYLE(LINE),THICKNESS(2)][])): >  display(Array([display(SP, CP, size=[400,400]),                display(newSP, CP, size=[400,400])]));         >

Download 2D_density_plot_with_Kernel_ac0.mw

I made no changes for efficiency to the earlier code. And I would assemble the replacement Matrix of plotting data differently -- using BlockCopy rather than an initializer with nested `if` calls. And one has to be very careful with the scaling of the short, thin lines relative to the plot view, or else weird visual effects can occur. This is not production level code. I'm really just trying to show that I like the idea of symbolsize=0.1 and so on. In these situations the plot may look less like it was drawn with crayons.

I don't see any information in your attachment that could be used to restore anything.

You wrote of some "given two dimensional functions" (emphasis mine).

If that means that you have those functions -- either as formulas or as means to compute at arbitrary x-y values -- then you might be able to sensibly use those to solve it as an optimization problem.

If you don't have those functions then note that the largest distance between the two represented curves may very well not occur at any of the known x[i],y[i] points. If the data points are fine enough then you might be satisfied by optimization of the difference of a pair or corresponding spline interpolations. You might be able to estimate error bounds on that (usually requiring some smoothness, etc).

Your Question lacks detail.

@mmcdara The link of your latest attachment isn't accessible for me, perhaps due to using "+" in the filename.

@sand15 If you wish to try smoothing via convolution then (for efficiency) you could look at the Convolution command from ImageTools.

I was also wondering a little about possibly using Threshold as a poor man's 2D bin separator. Of course it might be quite inefficient to perform multiple passes, in the case that one did not want just a few bins. But it's already there, and reasonably quick.

@mmcdara That link to R code showed some functionality involving contours overlaid with a density plot.

I've noticed in the past that neither 2D surfdata (working from 2D data, using dimension=2) nor densityplot (working from a bivariate function) allow for style=contour or style=surfacecontour, in a way similar to how they work for style=surface.

And it's not just the commands; the GUI renderer simply rejects STYLE(CONTOUR) as substructure in this case of a GRID & COLOR within a 2D PLOT structure.

But you can create an interpolant from the data, and generate the 2D contour plot from that.

Download 2D_density_plot_contour_ac.mw

I think that it would be an improvement if the GUI renderer allowed for this functionality directly, respecting STYLE(CONTOUR) or STYLE(PATCHCONTOUR).

Don't delete the full body of the original question. Here it is below:

Write a set of Maple procedures to perform the standard arithmetic operations on polynomials of a single variable with rational coefficients:

AddPolynomial(<poly1>,<poly2>) returns the sum SubPolynomial(<poly1>,<poly2>) returns the difference MultPolynomial(<poly1>,<poly2>) returns the product DivPolynomial(<poly1>,<poly2>) returns a list [<quotient>,<remainder>] GcdPolynomial(<poly1>,<poly2>), where deg(poly1) ≥ deg(poly2), returns the gcd

Begin by writing a utility procedure called CoeffList that converts a Maple polynomial into a list of its coefficients. For example, CoeffList(1/6*x^3 - 1/2*x^2 + x - 1) should return the list [1/6, -1/2, 1, -1].

Your arithmetic procedures should operate on, and return, lists that are in CoeffList format. Leading zero coefficients should be suppressed.

You can use Maple's facilities for arithmetic on rational numbers, but not for polynomial manipulation.

Following the same convention as Maple, GcdPolynomial should return a polynomial with integer coefficients if both of its arguments have integer coefficients, otherwise return a monic polynomial with rational coefficients.