Search code examples
pythonwolfram-mathematicanumerical-methodsdifferential-equations

defining a fuction to be used in solving equation from a data file


I am completely new to python and in fact any fundamental programming language, I use Mathematica for my all my symbolic and numeric calculations. I am learning to work with python and finding it really awesome! Here is a problem I am trying to solve but stuck without a clue! I have a data file for example

0.      1.    
0.01    0.9998000066665778    
0.02    0.9992001066609779    
...     ..
 

Which just the {t, Cos[2t]}. I want to define a function out of this data and use it in solving an equation in python. My Mathematica intuition tells me that I should define the function like:

    iFunc[x_] = Interpolation[iData, x]

and rest of the job is easy. for instance

NDSolve[{y''[x] + iFunc[x] y[x] == 0, y[0] == 1, y[1] == 0}, y, {x, 0, 1}]
    

Solves the equation easily. (I have not tried with more complicated cases though). Now how to do the job in python and also accuracy is an important issue for me. So, now I would like to ask two questions.

1. Is this the most accurate method in Mathematica?

2. And what is the equivalent of more accurate way to do the problem in python?

Here is my humble attempt to solve the problem (with a lot of input from StackOverflow) where the definition with cos(2t) works:

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

data = np.genfromtxt('cos2t.dat')

T = data[:,0]  #first column
phi = data[:,1]  #second column

f = interpolate.interp1d(T, phi)

tmin = 0.0# There should be a better way to define from the data
dt = 0.01
tmax = 2*np.pi
t = np.arange(tmin, tmax, dt) 

phinew = f(t)   # use interpolation function returned by `interp1d`

"""
def fun(z, t):
    x, y = z
    return np.array([y, -(cos(2*t))*x ])
"""    
def fun(z, t):
    x, y = z
    return np.array([y, -(phinew(t))*x ])    

sol1 = odeint(fun, [1, 0], t)[..., 0]
    
# for checking the plots        
plt.plot(t, sol1, label='sol')
plt.show()


        

*When I run the code with interpolated function from cos(2t) data, is not working...the error message tell

Traceback (most recent call last): File "testde.py", line 30, 
     in <module> sol1 = odeint(fun, [1, 0], t)[..., 0] 
File "/home/archimedes/anaconda3/lib/python3.6/site-packages/scip‌​y/integrate/odepack.‌​py", 
    line 215, in odeint ixpr, mxstep, mxhnil, mxordn, mxords) 
File "testde.py", 
    line 28, in fun return np.array([y, -(phinew(t))*x ])
TypeError: 'numpy.ndarray' object is not callable. 

I really can't decipher them. Please help...


Solution

  • To sub-question 2

    With

    t = np.arange(tmin, tmax, dt) 
    
    phinew = f(t)   # use interpolation function returned by `interp1d`
    

    equivalent to

    phinew = np.array([ f(s) for s in t])
    

    you construct phinew not as callable function but as array of values, closing the circle array to interpolation function to array. Use f which is a scalar function directly in the derivatives function,

    def fun(z, t):
        x, y = z
        return np.array([y, -f(t)*x ])