Plot Cassini ovals in Python using numpy and matplotlib

Question

I am trying to plot Cassini ovals in Python using these parametric equations for x,y.

https://mathcurve.com/courbes2d.gb/cassini/cassini.shtml

Here is my program.

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt


cnt = 20

plt.figure(figsize=(12,8))

for a in np.linspace(1,20,cnt):

    b = 0.9*a

    t = np.linspace(a + np.sqrt(a**2-b**2), a + np.sqrt(a**2+b**2),200)
    x = (b**4 - t**4)/(4*a*(t**2))
    y = np.sqrt(t**2 - ((x-a)**2))

    plt.plot(x, y)
    plt.plot(x, -y)

1. I am not sure if it works OK, I do get some nice plot but I am not sure if the Cassini ovals are full. Are they full/complete or just partially plotted? Btw, they just look like circles, not quite what I was hoping for.

2. When I increase cnt (say to 100), I start getting weird warnings RuntimeWarning: invalid value encountered in sqrt: y = np.sqrt(t**2 - ((x-a)**2)) I don't know why I am not getting these warnings when cnt is smaller. Why am I getting these for larger values of cnt?

How can I fix this code?

The first picture below is for cnt=20, the second one is for cnt=100.

The second one has these blank areas, I think that's where the quantity under the sqrt sign becomes negative (or something), and probably that's when I get these warnings. But why? I just programmed literally the two parametric equations. I am not quite sure what's wrong.

Answer 1

The ovals

Are they full/complete or just partially plotted?

I guess, they are full plotted.

Btw, they just look like circles, not quite what I was hoping for.

Your choice of parameters is b = 0.9a > 0, for which looking like circles is OK. I tried other relation between parameters, as in the link you provided, see example below.

I'm using Julia. So, Cassinian oval is

function cassinian(a, b)
    t = if a ? b
        range(a + sqrt(a^2 - b^2), a + sqrt(a^2 + b^2); length=200)
    else
        range(-a + sqrt(a^2 + b^2), a + sqrt(a^2 + b^2); length=200)
    end
    x = @. (b^4 - t^4) / (4 * a * t^2)
    y = @. sqrt(t^2 - (x-a)^2)
    return (x, y)
end

In your example b = 0.9a (a here goes from one to 20 with step 1, inclusively)

plt = plot(; title="b = 0.9a > 0", legend=nothing)
for a in 1:20
    x, y = cassinian(a, 0.9a)
    plot!(x, y)
    plot!(x, -y)
end
plt

If a ? b ? a?2 the "ovals" look no more like circles. Here a is one, b goes from 1 to ?2 (inclusively)

plt = plot(; title="a ? b ? a?2", legend=nothing)
for b in range(1, ?2; length=5)
    x, y = cassinian(1, b)
    plot!(x, y)
    plot!(x, -y)
end
plt

The errors

I also have errors for some pairs of a and b. For sure, they are caused by precision limits of floats. Even Julia's BigFloat can't handle this. More often, an error happens in calculation of first y, here a square root of a small negative value happens (about machine precision 1e-16). I've tried expansion of powers and this workaround for sqrt in generation of t vector, but that did not much help.

P.S.

Lovely curves <3. Did not hear about them, thank you!