** Functions on arrays of values
:PROPERTIES:
:date: 2013/02/27 14:49:49
:updated: 2013/03/06 19:38:28
:categories: python
:END:
It is common to evaluate a function for a range of values. Let us consider the value of the function $f(x) = \cos(x)$ over the range of $0 < x < \pi$. We cannot consider every value in that range, but we can consider say 10 points in the range. The func:numpy.linspace conveniently creates an array of values.
#+BEGIN_SRC python
import numpy as np
print np.linspace(0, np.pi, 10)
#+END_SRC
#+RESULTS:
: [ 0. 0.34906585 0.6981317 1.04719755 1.3962634 1.74532925
: 2.0943951 2.44346095 2.7925268 3.14159265]
The main point of using the mod:numpy functions is that they work element-wise on elements of an array. In this example, we compute the $\cos(x)$ for each element of $x$.
#+BEGIN_SRC python
import numpy as np
x = np.linspace(0, np.pi, 10)
print np.cos(x)
#+END_SRC
#+RESULTS:
: [ 1. 0.93969262 0.76604444 0.5 0.17364818 -0.17364818
: -0.5 -0.76604444 -0.93969262 -1. ]
You can already see from this output that there is a root to the equation $\cos(x) = 0$, because there is a change in sign in the output. This is not a very convenient way to view the results; a graph would be better. We use mod:matplotlib to make figures. Here is an example.
#+BEGIN_SRC python
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(0, np.pi, 10)
plt.plot(x, np.cos(x))
plt.xlabel('x')
plt.ylabel('cos(x)')
plt.savefig('images/plot-cos.png')
#+END_SRC
#+RESULTS:
[[./images/plot-cos.png]]
This figure illustrates graphically what the numbers above show. The function crosses zero at approximately $x = 1.5$. To get a more precise value, we must actually solve the function numerically. We use the function func:scipy.optimize.fsolve to do that. More precisely, we want to solve the equation $f(x) = \cos(x) = 0$. We create a function that defines that equation, and then use func:scipy.optimize.fsolve to solve it.
#+BEGIN_SRC python
from scipy.optimize import fsolve
import numpy as np
def f(x):
return np.cos(x)
sol, = fsolve(f, x0=1.5) # the comma after sol makes it return a float
print sol
print np.pi / 2
#+END_SRC
#+RESULTS:
: 1.57079632679
: 1.57079632679
We know the solution is \pi/2.