Solving an ordinary differential equation on a fixed grid (preferably in python)
<p>I have a differential equation of the form </p>

<p><code>dy(x)/dx = f(y,x)</code></p>

<p>that I would like to solve for <code>y</code>. </p>

<p>I have an array <code>xs</code> containing all of the values of <code>x</code> for which I need <code>ys</code>. </p>

<p>For only those values of <code>x</code>, I can evaluate <code>f(y,x)</code> for any <code>y</code>. </p>

<p>How can I solve for <code>ys</code>, preferably in python?</p>

<h2>MWE</h2>

<pre><code>import numpy as np

# these are the only x values that are legal
xs = np.array([0.15, 0.383, 0.99, 1.0001])

# some made up function --- I don't actually have an analytic form like this
def f(y, x):
    if not np.any(np.isclose(x, xs)):
        return np.nan
    return np.sin(y + x**2) 

# now I want to know which array of ys satisfies dy(x)/dx = f(y,x)
</code></pre>

<p>I think you should first solve ODE on a regular grid, and then interpolate solution on your fixed grid. The approximate code for your problem</p>

<pre><code>import numpy as np
from scipy.integrate import odeint
from scipy import interpolate

xs = np.array([0.15, 0.383, 0.99, 1.0001])

# dy/dx = f(x,y)
def dy_dx(y, x):
    return np.sin(y + x ** 2)

y0 = 0.0 # init condition
x = np.linspace(0, 10, 200)# here you can control an accuracy
sol = odeint(dy_dx, y0, x)
f = interpolate.interp1d(x, np.ravel(sol))
ys = f(xs)
</code></pre>

<p>But dy_dx(y, x) should always return something reasonable (not np.none).
Here is the drawing for this case
<a href="http://i.stack.imgur.com/j5RdJ.png" rel="nofollow"><img src="http://i.stack.imgur.com/j5RdJ.png" alt="enter image description here"></a></p>

python math mathematical-optimization numerical-methods differential-equations