Troubles with solving differental equations in Python

Question

I am trying to complete this, where I have to solve five ordinary differential equations using odeint and reproduce figures given in that task.

Here is my code:

import scipy as sp
import scipy.interpolate as ip
import numpy as np
import matplotlib.pyplot as pl

d = 8.64
Mu1 = 4.95*10**2
Mu2 = 4.95*10**(-2)
vs = 0.12
vd = 1.23
w = 10**(-3)
k1 = 2.19*10**(-4)
k2 = 6.12*10**(-5)
k3 = 0.997148
k4 = 6.79*10**(-2)

p0 = 1.00
sigmas0 = 2.01
sigmad0 = 2.23
alphas0 = 2.20
alphad0 = 2.26

hs = (sigmas0-(sigmas0**(2)-k3*alphas0*(2*sigmas0-alphas0))**(1/2))/k3
cs = (alphas0-hs)/2
ps = k4*(hs**2)/cs

t_points = [ 1000, 1850, 1950, 1980, 2000, 2050, 2080, 2100, 2120, 2150, 2225, 2300, 2500, 5000 ]
y_points = [ 0.0, 0.0, 1.0, 4.0, 5.0, 8.0, 10.0, 10.5, 10.0, 8.0, 3.5, 2.0, 0.0, 0.0 ]

t1 = np.array(t_points)
y1 = np.array(y_points)

new_length = 1000
new_t = np.linspace(t1.min(), t1.max(), new_length)
new_y2 = ip.pchip_interpolate(t1, y1, new_t)

pl.plot(t_points,y_points,'o', new_t,new_y2)
pl.show()

ft = sp.interpolate.interp1d(new_t, new_y2)

def equations(x, t1):

    p = x[0]
    alphad = x[1]
    alphas = x[2]
    sigmad = x[3] 
    sigmas = x[4]

    dpdt = (ps-p)/d + ft/Mu1
    dalphaddt = (1/vd)*(k2-w*(alphad-alphas))
    dalphasdt = (1/vs)*(w*(alphad-alphas)-k2)
    dsigmaddt = (1/vd)*(k1-w*(sigmad-sigmas))
    dsigmasdt = (1/vs)*(w*(sigmad-sigmas)-k1-(ps-p)/d*Mu2)

    return [dpdt, dalphaddt, dalphasdt, dsigmaddt, dsigmasdt]

solve =  sp.integrate.odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], t1)

It seems like this part:

dpdt = (ps-p)/d + ft/Mi1

is wrong and I have no idea how to solve it.

The error says:

TypeError: unsupported operand type(s) for /: 'interp1d' and 'float'.

Any ideas and advices are much appreciated.

EDIT: When I apply changes suggested by meowgoesthedog, I get error:

Traceback (most recent call last):

  File "<ipython-input-5-324757833872>", line 1, in <module>
    runfile('E:/Data/Project 2/project2.py', wdir='E:/Data/Project 2')

  File "D:\Programs\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile
    execfile(filename, namespace)

  File "D:\Programs\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "E:/Data/Project 2/project2.py", line 59, in <module>
    solve =  odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], t1)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\integrate\odepack.py", line 233, in odeint
    int(bool(tfirst)))

  File "E:/Data/Project 2/project2.py", line 51, in equations
    dpdt = (ps-p)/d + ft(t1)/Mu1

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\polyint.py", line 79, in __call__
    y = self._evaluate(x)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 664, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 696, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "

ValueError: A value in x_new is above the interpolation range.

`

Answer 1

According to interp1d's documentation:

ynew = f(xnew) # use interpolation function returned by interp1d

It returns a function / callable object which takes a value x and returns the interpolated value of f(x). In your case "x" = t:

dpdt = (ps-p)/d + ft(t1)/Mu1   # pass t1 to ft to obtain interpolated value

---

UPDATE

1. This new error is due to odeint sampling the function f(t) at values of t beyond the last value of t_points. This is necessary for error correction and there is no option to prevent odeint from doing so. However, we can instead extrapolate f(t) beyond the supplied samples, using InterpolatedUnivariateSpline:

   from scipy.interpolate import InterpolatedUnivariateSpline
   
   ...    
   
   ft = InterpolatedUnivariateSpline(t1, y1, k=1)
   

   As with interp1d, this returns a function with the same signature. However, after applying this fix the result becomes:

   

   Which is of course incorrect.

2. You have declared hs, cs, ps outside of the function as constants. In-fact they are functions of the alpha* and sigma* variables, so have to be evaluated during each call to equation:

   def equations(x, t):
   
       p = x[0]
       alphad = x[1]
       alphas = x[2]
       sigmad = x[3]
       sigmas = x[4]
   
       hs = (sigmas-(sigmas**(2)-k3*alphas*(2*sigmas-alphas))**(1/2))/k3
       cs = (alphas-hs)/2
       ps = k4*(hs**2)/cs
   
       dpdt = (ps-p)/d + ft(t)/Mu1
       dalphaddt = (1/vd)*(k2-w*(alphad-alphas))
       dalphasdt = (1/vs)*(w*(alphad-alphas)-k2)
       dsigmaddt = (1/vd)*(k1-w*(sigmad-sigmas))
       dsigmasdt = (1/vs)*(w*(sigmad-sigmas)-k1-(ps-p)/d*Mu2)
   
       return [dpdt, dalphaddt, dalphasdt, dsigmaddt, dsigmasdt]
   

   The result now matches the graph in the exercise... almost.

   

3. You passed t1 as the horizontal axis variable to odeint. It only has 14 elements which is too few for a smooth output. Pass new_t instead:

   solve = ig.odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], new_t)
   

   The result now exactly matches the expected one!