Approximating integral with infinite bounds

[saad khalid]

I'm trying to compute/estimate a rather complicated looking integral. Here is the code I'm trying to run, with the necessary constants defined:

var('t1,t2,u,w,k')
T = 1
m = 100
E = 1
v = 0
y=1
O = 1
integral(integral(integral(
   integral(integral(
     e^(-t1^2/T^2)*e^(-t2^2/T^2)*e^(I*O*t1)*
     e^(-I*O*t2)*e^(-I*E*y^2*(1 - v)*t1^2/2)*
      e^(-I*E*y^2*(1 - v)*t2^2/2)*e^(-I*k*y*(1 - u)*t1)*
      e^(I*k*y*(1 - v)*t2)*
      e^((1 + I)*(sqrt(E)*y*w*t1 + w*k/sqrt(E)))*
      e^((1 - I)*(sqrt(E)*y*u*t2 + u*k/sqrt(E)))*
      e^(-w^2/2)*e^(-u^2/2)*w^(-1/2 + I*m^2/(2*E))*
      u^(-1/2 - I*m^2/(2*E)), (u, 0, Infinity)), (w, 0,
     Infinity)), (t2, -Infinity, Infinity)), (t1, -Infinity,
   Infinity)), (k, -Infinity, Infinity))

I haven't been able to get a result from this code, it seems to run forever. I was hoping to be able to estimate the integral with some numerical methods, however I was having trouble getting a numerical integral set up properly. My first question is, can someone help me set up multivariable numerical integrals properly. I was trying something like

numerical_integral(x*y,(x,0,1),(y,0,1))

or

numerical_integral(numerical_integral(x*y,(x,0,1)),(y,0,1))

but neither seem to be the correct format, as they both give errors.

My second question is, can anyone give some advice on how to approximate such an integral, where it's multivariable and the bounds are at infinity? I don't really know where to start.

Thanks!

[saad khalid]

Also, I tried running the integral with finite bounds, and I get a giac error. Here is the code:

var('t1,t2,u,w,k')
T = 1
m = 100
E = 1
v = 0
y=1
O = 1
integral(integral(integral(
   integral(integral(
     e^(-t1^2/T^2)*e^(-t2^2/T^2)*e^(I*O*t1)*
     e^(-I*O*t2)*e^(-I*E*y^2*(1 - v)*t1^2/2)*
      e^(-I*E*y^2*(1 - v)*t2^2/2)*e^(-I*k*y*(1 - u)*t1)*
      e^(I*k*y*(1 - v)*t2)*
      e^((1 + I)*(sqrt(E)*y*w*t1 + w*k/sqrt(E)))*
      e^((1 - I)*(sqrt(E)*y*u*t2 + u*k/sqrt(E)))*
      e^(-w^2/2)*e^(-u^2/2)*w^(-1/2 + I*m^2/(2*E))*

      u^(-1/2 - I*m^2/(2*E)), (u, 0, 10)), (w, 0,
     10)), (t2, -10, 10)), (t1, -10,
   10)), (k, -10, 10))

Here is the error

RuntimeError: An error occurred running a Giac command:
INPUT:
sage20
OUTPUT:
:1: syntax error  line 1 col 31 at " in sage20:=int(sage16,sage17"Done",sage18w,sage190):;
 :1: syntax error  line 1 col 31 at " in sage20:=int(sage16,sage17"Done",sage18w,sage190):;
 "Done"

[Montgomery-Smith, Stephen]

Noting that there are lots of exp(-x^2) in the integrand, I would perform the numerical approximation using Monte-Carlo techniques with Gaussian pseudu-random numbers.

---

From: sage-s...@googlegroups.com <sage-s...@googlegroups.com> on behalf of saad khalid <saad...@gmail.com>
Sent: Wednesday, February 26, 2020 9:41:19 PM
To: sage-support
Subject: [sage-support] Re: Approximating integral with infinite bounds

 

--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/2ee464b7-bb34-4365-8c5e-268cf5d655fb%40googlegroups.com.

[saad khalid]

Hi, I was trying to use the monte-carlo integration in Sage. However, I couldn't get it to work because there are complex numbers involved. Whenever I would try to integrate something like e^(i*x), the monte-carlo integrator would complain about it not being an integral of real numbers. Do you have any suggestion or reference on how to approach this?

Thanks!

To unsubscribe from this group and stop receiving emails from it, send an email to sage-s...@googlegroups.com.

[Montgomery-Smith, Stephen]

Split it into it's real and imaginary parts maybe?

---

From: sage-s...@googlegroups.com <sage-s...@googlegroups.com> on behalf of saad khalid <saad...@gmail.com>

Sent: Thursday, February 27, 2020 2:51:07 PM
To: sage-support
Subject: Re: [sage-support] Re: Approximating integral with infinite bounds

 

To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/02dc46ce-545a-4ad0-849b-e8668ec6282a%40googlegroups.com.

[Emmanuel Charpentier]

This question would have been more properly posed on ask.sagemath.org...

The firs argument of numerical integral must be a function of one (real) variable returning a real. It is indeed possible to nest such integrals, but the second one nests the first one, hence long computation times. Examples:

sage: (e^-(x^2/2)).integrate(x,-oo,oo)

sqrt(2)*sqrt(pi)

sage: (e^-(x^2/2)).integrate(x,-oo,oo).n()

2.50662827463100

sage: %time numerical_integral(lambda x1:e^-(x1^2/2),-oo,oo)

CPU times: user 29.8 ms, sys: 21 µs, total: 29.9 ms

Wall time: 29.8 ms

(2.5066282746309994, 2.5512942303338327e-08)

Both precision and computation time are acceptable. But consider:

sage: var("y")

y

sage: (e^-((x^2+y^2)/2)).integrate(y,-oo,oo).integrate(x,-oo,oo)

2*pi

sage: (e^-((x^2+y^2)/2)).integrate(y,-oo,oo).integrate(x,-oo,oo).n()

6.28318530717959

sage: %time numerical_integral(lambda x2:numerical_integral(lambda x1:e^-((x1^2+x

....: 2^2)/2),-oo,oo)[0],-oo,oo)

CPU times: user 4.05 s, sys: 32 ms, total: 4.08 s

Wall time: 4.08 s

(6.283185307182902, 6.39626298408729e-08)

Acceptable precision, coputation time questionable. And this is only for a double integral of a "well behaved" real function.

Furthermore, as pointed out, the values of your function are complex. You must, at each step, separate real and imaginary part...

Therefore, I doubt that this solution is viable.

Monte-Carlo integration has been suggested. This solution has been well-explored in some R packages dedicated to Bayesian statistics, but the computation of the =value of the integral is still a delicate question (and, being only of secondary interest to applied statistics, not as well  explored).

You may try to separate the real and imaginary part of your function, and sample from these part. Look at Stan and the companion R package bridgesampling (which could compute the integral from a pseudo-random sample obtained via Stan).

HTH,

[Vincent Delecroix]

Le 29/02/2020 à 21:41, Emmanuel Charpentier a écrit :
>
> You may try to separate the real and imaginary part of your function, and

> sample from these part. Look at Stan<https://mc-stan.org/> and the
> companion R package bridgesampling
> <https://cran.r-project.org/web/packages/bridgesampling/index.html> (which

> could compute the integral from a pseudo-random sample obtained via Stan).

You can do Monte-Carlo in Sage (via GSL)

sage: f(x,y,z,t) = x^2*y*t + 3*x*y^2 - x*y*z + t
sage: monte_carlo_integral(f, [-1,-1,-1,-1], [1,1,1,1], 10^6)
(0.0155081492100384, 0.015914524378217675)

But you can not specify manually the order of the variables
which is a bit confusing.

Vincent

[Simon King]

Hi Emmanuel,

On 2020-02-29, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:
> This question would have been more properly posed on ask.sagemath.org...

Why?

I for one totally dislike ask.sagemath.org and would never post a
question or an answer there. It is of course a matter of taste. But it
is certainly not appropriate to discourage people using sage-support,
IMHO.

Best regards,
Simon

[Fernando Gouvea]

Some years ago in a book review, David Roberts had the idea of plotting an algebraic curve using the transformation  (u,v) = (x,y)/(r² + x² + y²)^1/2, which transforms the plane into a circle and makes it easy to visualize the projective completion of the curve. You can see some of his plots at https://www.maa.org/press/maa-reviews/rational-algebraic-curves-a-computer-algebra-approach

I’d love to do this kind of plot for my students. Can anyone offer help on how to do it with Sage? (Of course the dream scenario would be to add this option to the plot method for curves...) 

I’ve been using implicit_plot for most of my examples, which seems to be equivalent of using C.plot() when C is a curve.

Thanks,

Fernando

--

==================================================================
Fernando Q. Gouvea                                          Editor, MAA Reviews
Dept of Mathematics and Statistics                     http://www.colby.edu/~fqgouvea
Colby College                                                    http://www.maa.org/press/maa-reviews
Mayflower Hill 5836              
Waterville, ME 04901         

A training in mathematics is a prerequisite today for work in almost
any scientific field, but even for those who are not going to become
scientists, it is essential because, if it is only through speech that
we can understand what freedom means, only through mathematics
can we understand what necessity means.
  -- W. H. Auden

[Emmanuel Charpentier]

Le samedi 29 février 2020 23:06:36 UTC+1, Simon King a écrit :

Hi Emmanuel,

On 2020-02-29, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:
> This question would have been more properly posed on ask.sagemath.org...

Why?

As I understand our current documentation :

sage-support is oriented towards reporting problems in Sage itself (bugs, problematic results, unimplmented or ill-implemented features, etc...). You go there when you (mathematically) know what you want to do but Sage doesn't do it  as (you think) it should.

ask.sagemath.org is more oriented towards mathematical questions: "how can I solve this particular problem with Sage ?".

 

I for one totally dislike ask.sagemath.org

I am not exceedingly fond of its infantilizing user interface either...

 

and would never post a  question or an answer there.

People (often Sage newcomers) do ask questions here (probably because they are told to by our documentation...). They may deserve an answer (even if this answer is "We won't do your homework for you", as it happens too often...).

 

It is of course a matter of taste. But it is certainly not appropriate to discourage people using sage-support, IMHO.

Could you propose an update to our documentation ? 

Best regards,
Simon

[Fernando Gouvea]

(Repost, as the first message seems to have ended up in the wrong thread.)

Some years ago in a book review, David Roberts had the idea of plotting an algebraic curve using the transformation  (u,v) = (x,y)/(r² + x² + y²)^1/2, which transforms the plane into a circle and makes it easy to visualize the projective completion of the curve. You can see some of his plots at https://www.maa.org/press/maa-reviews/rational-algebraic-curves-a-computer-algebra-approach

I’d love to do this kind of plot for my students. Can anyone offer help on how to do it with Sage? (Of course the dream scenario would be to add this option to the plot method for curves...) 

I’ve been using implicit_plot for most of my examples, which seems to be equivalent of using C.plot() when C is a curve.

Thanks,

Fernando

-- 
=============================================================
Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea
Carter Professor of Mathematics
Dept. of Mathematics and Statistics
Colby College              
5836 Mayflower Hill        
Waterville, ME 04901       

Why is the alphabet in that order? Is it because of that song?

[David Joyner]

On Sat, Feb 29, 2020 at 5:29 PM Fernando Gouvea <fqgo...@colby.edu> wrote:

Some years ago in a book review, David Roberts had the idea of plotting an algebraic curve using the transformation  (u,v) = (x,y)/(r² + x² + y²)^1/2, which transforms the plane into a circle and makes it easy to visualize the projective completion of the curve.

 

You can see some of his plots at https://www.maa.org/press/maa-reviews/rational-algebraic-curves-a-computer-algebra-approach

These look like 2d parametric plots.

https://doc.sagemath.org/html/en/reference/plotting/sage/plot/plot.html#sage.plot.plot.parametric_plot

To combine two (or more) par. plots on the same axis, just add the corresponding graphics objects.

There are examples on the pages linked to above.

 

I’d love to do this kind of plot for my students. Can anyone offer help on how to do it with Sage? (Of course the dream scenario would be to add this option to the plot method for curves...) 

I’ve been using implicit_plot for most of my examples, which seems to be equivalent of using C.plot() when C is a curve.

Thanks,

Fernando

--

==================================================================
Fernando Q. Gouvea                                          Editor, MAA Reviews
Dept of Mathematics and Statistics                     http://www.colby.edu/~fqgouvea
Colby College                                                    http://www.maa.org/press/maa-reviews
Mayflower Hill 5836              
Waterville, ME 04901         

A training in mathematics is a prerequisite today for work in almost
any scientific field, but even for those who are not going to become
scientists, it is essential because, if it is only through speech that
we can understand what freedom means, only through mathematics
can we understand what necessity means.
  -- W. H. Auden

--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/CAPON7s8KPbT-cLo_4FHJN-WqoqHi82%3DmgS1HBAD8SJqDxqkFDg%40mail.gmail.com.

[Fernando Gouvea]

Parametric plots won't work for general algebraic curves.

But I'm also not sure how to implement the transformation into a circle. I'll look at the documentation for plots.

Fernando

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/CAEQuuAUukezRFvtNBmhb6VPfkT0AjvbbGE55%2BM7xj90ZhLk2Lw%40mail.gmail.com.

-- 
=============================================================
Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea

Carter Professor of Mathematics
Dept. of Mathematics and Statistics
Colby College              
5836 Mayflower Hill        
Waterville, ME 04901       

There is nothing mysterious, as some have tried to maintain, about the
applicability of mathematics. What we get by abstraction from
something can be returned.
  -- R. L. Wilder, Introduction to the Foundations of Mathematics.

[Fernando Gouvea]

Here's what I ended up trying, with r=3:

var('x y u v')
x=u*sqrt(9/(1-u^2-v^2))
y=v*sqrt(9/(1-u^2-v^2))
implicit_plot(y^2-x^3+x==0,(u,-1,1),(v,-1,1))

That gives an error:

/opt/sagemath-8.9/local/lib/python2.7/site-packages/sage/ext/interpreters/wrapper_rdf.pyx in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.__call__ (build/cythonized/sage/ext/interpreters/wrapper_rdf.c:2237)()
     74         for i from 0 <= i < len(args):
     75             self._args[i] = args[i]
---> 76         return self._domain(interp_rdf(c_args
     77             , self._constants
     78             , self._py_constants

ValueError: negative number to a fractional power not real

Is there some way to tell implicit_plot to stay inside u^2+v^2\leq 1? Or to ignore complex values?

The equivalent code seems to give the correct graph in Mathematica.

Fernando

-- 
=============================================================
Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea

Carter Professor of Mathematics
Dept. of Mathematics and Statistics
Colby College              
5836 Mayflower Hill        
Waterville, ME 04901       

If little else, the brain is an educational toy.
  -- Tom Robbins

[Dima Pasechnik]

On Tue, Mar 3, 2020 at 8:20 PM Fernando Gouvea <fqgo...@colby.edu> wrote:
>
> Here's what I ended up trying, with r=3:
>
> var('x y u v')
> x=u*sqrt(9/(1-u^2-v^2))
> y=v*sqrt(9/(1-u^2-v^2))
> implicit_plot(y^2-x^3+x==0,(u,-1,1),(v,-1,1))
>
> That gives an error:
>
> /opt/sagemath-8.9/local/lib/python2.7/site-packages/sage/ext/interpreters/wrapper_rdf.pyx in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.__call__ (build/cythonized/sage/ext/interpreters/wrapper_rdf.c:2237)()
> 74 for i from 0 <= i < len(args):
> 75 self._args[i] = args[i]
> ---> 76 return self._domain(interp_rdf(c_args
> 77 , self._constants
> 78 , self._py_constants
>
> ValueError: negative number to a fractional power not real
>
> Is there some way to tell implicit_plot to stay inside u^2+v^2\leq 1? Or to ignore complex values?

I'd just change the limits of u and v to make the rectangle of the
values you plot in, anyway,
to well stay inside the unit circle.

> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/b1b0c01c-7dfd-0e17-a3d4-61012ab66d8b%40colby.edu.

[Dima Pasechnik]

On Tue, Mar 3, 2020 at 9:15 PM Dima Pasechnik <dim...@gmail.com> wrote:
>
> On Tue, Mar 3, 2020 at 8:20 PM Fernando Gouvea <fqgo...@colby.edu> wrote:
> >
> > Here's what I ended up trying, with r=3:
> >
> > var('x y u v')
> > x=u*sqrt(9/(1-u^2-v^2))
> > y=v*sqrt(9/(1-u^2-v^2))
> > implicit_plot(y^2-x^3+x==0,(u,-1,1),(v,-1,1))
> >
> > That gives an error:
> >
> > /opt/sagemath-8.9/local/lib/python2.7/site-packages/sage/ext/interpreters/wrapper_rdf.pyx in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.__call__ (build/cythonized/sage/ext/interpreters/wrapper_rdf.c:2237)()
> > 74 for i from 0 <= i < len(args):
> > 75 self._args[i] = args[i]
> > ---> 76 return self._domain(interp_rdf(c_args
> > 77 , self._constants
> > 78 , self._py_constants
> >
> > ValueError: negative number to a fractional power not real
> >
> > Is there some way to tell implicit_plot to stay inside u^2+v^2\leq 1? Or to ignore complex values?
>
> I'd just change the limits of u and v to make the rectangle of the
> values you plot in, anyway,
> to well stay inside the unit circle.

A better approach would be to use a custom scale --- note that
implicit_plot supports "scale" option,
which, I suppose, gets passed to matplotlib in some way for plotting,
I just don't know the details enough to
see how it can be implemented.
https://matplotlib.org/3.1.1/gallery/scales/custom_scale.html#sphx-glr-gallery-scales-custom-scale-py

[Fernando Gouvea]

The whole point of this is to show the behavior of the curve near infinity, so changing the limits is not an option.

Fernando

[Augustin Lefèvre]

A caveat is that at the boundary, the mapping you describe becomes non differentiable (the determinant of the differential blows up to infinity),

so it's going to be painful for implicit_plot to work.

That being said, the following tweak runs ok but it's not exactly what you describe.

var('x y u v')

x=u*sqrt(9/abs(1-u^2-v^2))
y=v*sqrt(9/abs(1-u^2-v^2))
pl=implicit_plot(y^2-x^3+x==0,(u,-1,1),(v,-1,1))
pl+=implicit_plot(u^2+v^2==1,(u,-1,1),(v,-1,1))
pl.show()



This other tweak raises an error, I don't see why :

var('x y u v')

x=u*sqrt(9/max(1-u^2-v^2,0))
y=v*sqrt(9/max(1-u^2-v^2,0))
pl=implicit_plot(y^2-x^3+x==0,(u,-1,1),(v,-1,1))
pl+=implicit_plot(u^2+v^2==1,(u,-1,1),(v,-1,1))
pl.show()

>> ValueError: negative number to a fractional power not real

[Dima Pasechnik]

On Tue, Mar 3, 2020 at 10:10 PM Fernando Gouvea <fqgo...@colby.edu> wrote:
>
> The whole point of this is to show the behavior of the curve near infinity, so changing the limits is not an option.

just paste together a number of rectangles where (u,v) stay inside the
unit circle.
(yes, this would need writing a loop, ideally)

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/6117cc36-6771-179c-2887-38c05305fd58%40colby.edu.

[Dima Pasechnik]

even better:

sage: var('x y u v r phi')
....: u=r*cos(phi)
....: v=r*sin(phi)
....: x=u*sqrt(9/(1-r^2))
....: y=v*sqrt(9/(1-r^2))
....: implicit_plot(y^2-x^3+x==0,(r,0,999/1000),(phi,-pi,pi))

[Fernando Gouvea]

Nice idea. Thanks.

Fernando

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/CAAWYfq1gnU%2BH55XmFEueQvPe8JTaiRSJ7G9k0cc8hSbdY2S9Zw%40mail.gmail.com.

[Fernando Gouvea]

But no, it doesn't work, since it gives a rectangular plot instead of one in polar coordinates. But maybe we are closer.

I still think implicit_plot should be smarter about values that do not make sense.

Fernando

[Dima Pasechnik]

On Wed, Mar 4, 2020 at 12:20 AM Fernando Gouvea <fqgo...@colby.edu> wrote:
>
> But no, it doesn't work, since it gives a rectangular plot instead of one in polar coordinates. But maybe we are closer.

I looked at the labels on the axes, and they do match the ranges of r
and phi, so I don't udnerstand
how it's possible.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/f7ea288f-80fc-478f-a89c-f514612b6802%40colby.edu.

[Dima Pasechnik]

The easiest way is to use Python functions rather than symbolic ones;
define a function that is 1 outside the unit disk, and implicitly plot it.

sage: def f_uv(u,v):
....:     if u^2+v^2>=1:
....:         return 1
....:     else:
....:         x=u*sqrt(9/(1-u^2-v^2))
....:         y=v*sqrt(9/(1-u^2-v^2))
....:         return y^2-x^3+x
....: implicit_plot(f_uv,(u,-1,1),(v,-1,1))

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/b1b0c01c-7dfd-0e17-a3d4-61012ab66d8b%40colby.edu.
>
> --
> =============================================================
> Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea
> Carter Professor of Mathematics
> Dept. of Mathematics and Statistics
> Colby College
> 5836 Mayflower Hill
> Waterville, ME 04901
>
> If little else, the brain is an educational toy.
>   -- Tom Robbins
>
> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/6117cc36-6771-179c-2887-38c05305fd58%40colby.edu.
>
> --
> =============================================================
> Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea
> Carter Professor of Mathematics
> Dept. of Mathematics and Statistics
> Colby College
> 5836 Mayflower Hill
> Waterville, ME 04901
>
> If little else, the brain is an educational toy.
>   -- Tom Robbins
>
> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.

[Fernando Gouvea]

This works, in the sense that there's no error. One does get a bunch of extraneous points near the boundary of the disk. It's as if plot_points were trying to connect the point at (0,1) and the point at (0,-1) along the circle, even though f_uv is 1 on the circle.

Strangely, they occur only on the right hand side (i.e., positive u, not negative u). I tried setting plot_points to be 500, but the bad points don't go away. Changing the curve to y^2-x^3+x-1=0 doesn't make them go away either.

Fernando

On 3/5/2020 8:22 AM, Dima Pasechnik wrote:

The easiest way is to use Python functions rather than symbolic ones;
define a function that is 1 outside the unit disk, and implicitly plot it.

sage: def f_uv(u,v):
....:     if u^2+v^2>=1:
....:         return 1
....:     else:
....:         x=u*sqrt(9/(1-u^2-v^2))
....:         y=v*sqrt(9/(1-u^2-v^2))
....:         return y^2-x^3+x
....: implicit_plot(f_uv,(u,-1,1),(v,-1,1))

>

> On Tue, Mar 3, 2020 at 8:20 PM Fernando Gouvea <fqgo...@colby.edu> wrote:

-- 
=============================================================
Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea
Carter Professor of Mathematics
Dept. of Mathematics and Statistics
Colby College              
5836 Mayflower Hill        
Waterville, ME 04901       

We now face a choice between Christ and nothing, because Christ has
claimed everything so that renouncing him can only be nihilism.
  -- Peter Leithart

[Dima Pasechnik]

On Thu, Mar 5, 2020 at 2:32 PM Fernando Gouvea <fqgo...@colby.edu> wrote:
>
> This works, in the sense that there's no error. One does get a bunch of extraneous points near the boundary of the disk. It's as if plot_points were trying to connect the point at (0,1) and the point at (0,-1) along the circle, even though f_uv is 1 on the circle.
>
> Strangely, they occur only on the right hand side (i.e., positive u, not negative u). I tried setting plot_points to be 500, but the bad points don't go away. Changing the curve to y^2-x^3+x-1=0 doesn't make them go away either.
>

the reason is that implicit_plot attempts to approximate the function
it assumes continuous, so if it's negative inside, but near, the
boundary, and positive nearby, but outside, then a fake zero is being
drawn very close to the boundary.

That's why it should be better to create a plot in polar coordinates
and then transform it.

> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/bc81b0c4-62ff-65b9-a7f2-e995f26f9234%40colby.edu.

[Dima Pasechnik]

In fact, substituting x and y directly into the equation of the curve
to plot, and clearing denominators,
produces something pretty good,IMHO:

implicit_plot(v^2*3*sqrt(1-u^2-v^2)-u^3*9+u*(1-u^2-v^2),(u,-1,1),(v,-1,1))

[Dima Pasechnik]

More conceptually, one can use, with care, Sage's substitution facilities:
sage: var('u v x y t');
sage: f=y^2-x^3+x
sage: fs=(f.subs(x=u*3*t^(-1/2),y=v*3*t^(-1/2))*t^(3/2)).expand() #
only works with extra variable t
sage: implicit_plot(fs.subs(t=1-u^2-v^2),(u,-1,1),(v,-1,1))

[Fernando Gouvea]

Yes, and I should have thought of that!

Fernando