``scipy.odeint`` giving different answer than analytical

Question

I was using scipy.integrate.odeint function , the ode is

$$\frac{y\ dx - x\ dy}{(x+y)^2} + dy = dx$$

with solution

$$y^2 - x^2 - y = c (x + y)\ .$$

Solving it via odeint

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
def dy_dx(y, x):
  return ((x+y)**2 - y)/((x+y)**2 - x)

xs = np.linspace(1,100,100)
y0 = 1
ys = odeint(dy_dx, y0, xs ,atol = 2)

ys = np.array(ys).flatten()

plt.plot(xs, ys)
plt.xlabel("x")
plt.ylabel("y")

which certainly doesn't match the answer. What could be going wrong?

Answer 1

Nothing is going wrong.

• Your differential equation is: $$\frac{\mathrm{d}y}{\mathrm{d}x} = \frac{(x+y)^2 - y}{(x+y)^2 - x}$$ Obviously if $x=y$, you have $\frac{\mathrm{d}y}{\mathrm{d}x} = 1$, i.e., $x$ grows exactly as “fast” as $y$. Thus, if $x=y$ for your initial conditions, it will stay that way.

• Your initial conditions are $x=y=1$ and your numerical result looks like $x=y$.

• For $c=-\tfrac{1}{2}$ and $c=\tfrac{1}{2}$, respectively, $x=y$ solves your analytical solutions (“solved answers”).