Search code examples
python-2.xodeodeint

odeint doesn't give me what expected

I'm starting with ODE and I have a python code that first tries to solve a DOE with the odeint function and then compares that solution with the solution for the ODE that you computed on your own (which you must insert into the code). However they don't quite coincide.

I've tried with software programs and it seems that the solution I computed manually is ok so I wonder why the odeint function doesn't give me what I expected.

import numpy as np
import matplotlib.pyplot as plt

# Define a function which calculates the derivative
def dy_dx(y,x):
    dydx = 2.0*x*(y-1)/(1-x**2) #
    return dydx

xs = np.linspace(-4,4,100) 

for i in range(-10,10,1):
    y0 = i
    ys = odeint(dy_dx, y0, xs) 
    ys = np.array(ys).flatten() 
    plt.plot(xs, ys);

# SECOND PART:

y_exact =  1+(y0)/(1-x**2) # manually computed solution

y_difference = ys - y_exact 
plt.subplot(2, 1, 1)
plt.plot(xs, ys, xs, y_exact, "--"); 
plt.title("Difference between exact solution and computed solution");

So I added the "range() part" to see how it varies with different initial conditions but all of them are along the axis x=-1

The solution I found manually has a limit there but it's not just a line, as you can see if you run the second part or visit https://www.symbolab.com/solver/ordinary-differential-equation-calculator/2x%5Cleft(y-1%5Cright)dx%2B%5Cleft(x%5E%7B2%7D-1%5Cright)dy%3D%200

I just wonder where the mistake is or why odeint gives that as result.difference between odeint and my solution

I should also add that the ODE might be a little bit weird because you get absolute values when integrating. That might be related.

Solution

  • Your initial condition for the numerical solution is y(-4)=y0, as odeint takes the first point of the time interval as initial time. Correspondingly you would have to change your exact solution to

    y_exact =  1+((y0-1)*(1-xs[0]**2))/(1-xs**2)

    as you can check with the tool of your choice, for instance WA y'(x)=2*x*(y(x)-1)/(1-x^2), y(a)=b.

    As you also observed, your ODE is singular at x=-1 and x=1, so that any solution only has as domain any of the 3 sub-intervals created by these points, namely the one containing the initial time. Also, numerical methods do not work nicely close to singularities, so you would have to restrict your integration domain to [-4, -1-delta] for some small but not too small delta. Or if you really wanted to explore the middle interval with initial time at 0, you need to perform two integrations, one from 0 to -1+delta and one from 0 to 1-delta.

    If you consider the related differential equation without real singularities

    y'(x) = -2*x*(y-1)/(1+x^2),  y(x0) = y0

    with exact solution

    y(x) = 1 + (y0-1)*(1+x0^2)/(1+x^2)

    implemented via the modified code

    import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define a function which calculates the derivative
def dy_dx(y,x): return  -2.0*x*(y-1)/(1+x**2) #

# Flow function
def phi(x,x0,y0): return 1 + (y0-1)*(1+x0**2)/(1+x**2)

xs = np.linspace(-4,4,100) 

for i in range(-10,10,2):
    y0 = 1+0.1*i
    ys = odeint(dy_dx, y0, xs) 
    ys = ys.flatten() 
    plt.plot(xs, phi(xs,xs[0],y0),c='lightgray', lw=4);
    plt.plot(xs, ys, lw=1, label='$y_0=%.2f$'%y0);
plt.grid(); plt.xlim(-5.5,4.2); plt.legend(); plt.show()

    you get the following plot were you find the (colored) numerical solution nicely centered inside the (underlayed gray) exact solution. enter image description here