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 p

Solution 1:

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