CODEWITHSUNDEEP
pythonode
Joe
Joe

Reputation: 189

Solving an ode y'=f (x) with numerical values of f (x) but without analitical expresion

I want to solve an ODE numerically in python, like y'=f(x) (with boundary condition y(0)=0). I don't know what the analytical expression of function f(x), instead I have a set of points (data) of this function for the domain where I want to solve the differential equation.

I have tried with odeint. But this method works when you know the explicit analytical expression for f(x), which is not my case. In particular, I have the following code with the function f(x) in an array (For simplicity I consider a known f(x), but in my real problem this array f(x) comes from a numerical simulation without known analytical expression).

The following code doesn't work, since odeint thinks that I have y'=x, and not my values f(x):

from scipy.integrate import odeint
import numpy as np

def dy_dx(y, f):
    return f #it doesn't work!!


xs = np.linspace(0,10,100)

f = np.sin(xs)*np.exp(-0.1*xs) #data of the function f, but in my real problem I DON'T KNOW THE ANALITICAL SOLUTION! JUST ONLY the points

ys = odeint(dy_dx, 0.0, xs)

There must be something in Python that can solve this. Basically you are solving the ode numerically and you know what the values of f(x) in the domain of the ode.

Upvotes: 1

Views: 301

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You should be able to solve this using the quadrature routines of scipy.integrate. If you really want to use the complicated form, you have to use interpolation, for instance as in

from scipy.integrate import odeint
from scipy.interpolate import interp1d
import numpy as np

xs = np.linspace(0,10,100+1);
fs = np.sin(xs)*np.exp(-0.1*xs) # = Imag( exp((1j-0.1)*x) )
# the exact anti-derivative of f is 
# F = Imag( (exp((1j-0.1)*x)-1)/(1j-0.1) )
#   = Imag( -(1j+0.1)*(exp((1j-0.1)*x)-1)/(1.01) )
#   = 1/1.01 - exp(-0.1*x)/1.01 * ( cos(x) + 0.1*sin(x) )

def IntF(x): return (1-np.exp(-0.1*x)*(np.cos(x)+0.1*np.sin(x)))/1.01 

f = interp1d(xs, fs, kind="quadratic", fill_value="extrapolate")

def dy_dx(y, x):
    return f(x) 

ys = odeint(dy_dx, 0.0, xs)

for x,y in zip(xs, ys): print "%8.4f %20.15f %20.15f"%(x,y,IntF(x))

with the first 10 lines

    x          interpolated           exact
 --------------------------------------------------
  0.0000    0.000000000000000    0.000000000000000
  0.1000    0.004965420470493    0.004962659238991
  0.2000    0.019671988500299    0.019669801188631
  0.3000    0.043783730081358    0.043781529336000
  0.4000    0.076872788780423    0.076870713937278
  0.5000    0.118430993242631    0.118428986914274
  0.6000    0.167875357240100    0.167873429717074
  0.7000    0.224555718642310    0.224553873611032
  0.8000    0.287762489870417    0.287760727322230
  0.9000    0.356734939606963    0.356733243391002
  1.0000    0.430669760236151    0.430668131955269

Upvotes: 3

Related Questions