[Nonsense question stemming from a simple typo]

Question

My confusion was just based on a typo in my code - nothing interesting to learn here! I simply plotted the wrong function, see my comment to the accepted answer.

Here are some seemingly harmless Python code

import matplotlib.pyplot as plt 
import numpy as np 
  
delta = 0.02    
    
x = np.arange(-10, 10, delta) 
y = np.arange(-10, 10, delta) 
[X, Y] = np.meshgrid(x, y) 
  
fig, ax = plt.subplots() 
  
Z = X**2+2*X*Y+Y**2
# for comparison later use this:
#Z = 3.14*X**2+0.6*Y**2 
  
ax.contour(X, Y, Z)   
ax.set_xlabel('x') 
ax.set_ylabel('y') 
  
plt.show() 

and its output:

The plotted contour lines are wrong! I am doing the contour plot of a quadratic function and the contour lines should be ellipses. The shown lines are not just seemingly lines, which at a bigger scale would become ellipses. I have at least two reasons for saying so:

1. Here is the Wolfram Alpha contour plot of the same function at a similar scale. It shows the ellipsis that should be there.

2. I calculated by hand the corresponding diagonal quadratic form (according to the Principal Axis theorem): It is, with a bit of rounding, Z = 3.14*X**2+0.6*Y**2.

This a function which first rotates a given vector (X,Y) by a certain angle and then applies the function Z = X**2+2*X*Y+Y**2 used in the code. Thus the contour plots of the two functions should just differ by a rotation.

But inserting Z = 3.14*X**2+0.6*Y**2 in the above code, in place of the previous Z, produces a perfectly fine ellipsis:

My questions:

1. What is the explanation for his behaviour?
2. What can I do to get an accurate contour plot of the first function?

Answer 1

The matplotlib contour plot is correct.

A) We can look at the 3D projection:

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10,10))
ax = fig.gca(projection='3d')

delta = 2  
    
x = np.arange(-100, 100, delta) 
y = np.arange(-100, 100, delta) 
X, Y = np.meshgrid(x, y)   
  
Z = X**2 + Y**2 + 2*X*Y
ax.plot_surface(X, Y, Z)

plt.show()

Sample output:

Yup, lines parallel to x = -y.

B) Let's look at the mathematical representation for x+y=a, i.e., y=a-x

Then it follows:
z = x² + y² + 2xy
= x² + (a-x)² + 2x(a-x)
= x² + a² - 2ax + x² + 2ax - 2x²
= a²

And this is what we see in the surface plot, as well as in the contour plot - for lines parallel to x=-y, the z-value is independent of x (or y) and just depends on a (and is symmetrical for positive and negative a-values).

C) Obviously, the contour plots differ for a in x²+2xy+ay² (as you expected). But because we can, let's add this little example:

import matplotlib.pyplot as plt 
import numpy as np 
  
delta = 0.02  
    
x = np.arange(-10, 10, delta) 
y = np.arange(-10, 10, delta) 
X, Y = np.meshgrid(x, y) 
  
fig, axes = plt.subplots(3, 3, figsize=(10,10)) 
  
Z = X**2 + Y**2 + 3*Y*X
# for comparison later use this:
#Z = 3.14*X**2+0.6*Y**2 

for i, ax in enumerate(axes.flatten()):
    a=i-2
    Z = X**2 + a* Y**2 + 2*X*Y
    cs = ax.contour(X, Y, Z)
    ax.clabel(cs, inline=True, fontsize=10)
    ax.set_title("x²+2xy+ay² for a="+str(a))
  
plt.tight_layout()
plt.show() 

Sample output: