Numerical python book practice

May 16, 2016

In [1]: %lsmagic
Out[1]: Available line magics:
%alias %alias_magic %autocall
Available cell magics:
%%! %%HTML %%SVG %%bash

%automagic

%%capture

%autosave

%%debug

%%file

Automagic is ON, % prefix IS NOT needed for line magics.

In [2]: %run test.py
ok! you have succesfully run this script

In [3]: %reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? n
Nothing done.

In [4]: %pylab notebook

Populating the interactive namespace from numpy and matplotlib

In [5]: data=array([[1,2],
[3,4],
[5,6]])
In [6]: data.shape
Out[6]: (3, 2)
In [7]: X=Y=linspace(1,10)
In [8]: type(X)
Out[8]: numpy.ndarray
1

%bookmark

%%html

%cat

%c

%%javascrip

In [9]: x,y=meshgrid(X,Y)
In [10]: from sympy import *
x=symbols('x')
func=lambda x:x**2 +2*x -3
func(x)
Out[10]: x**2 + 2*x - 3
In [11]: %reset
Once deleted, variables cannot be recovered. Proceed (y/[n])? n
Nothing done.

In [12]: %pylab notebook

from sympy import *
init_printing()
Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['test', 'floor', 'beta', 'sqr
`%matplotlib` prevents importing * from pylab and numpy
In [13]: x,y,z=symbols('x,y,z')
z=(x+y)**2
z
Out[13]:
(x + y)2
In [14]: x=y=linspace(1,3)
X,Y=meshgrid(x,y)
Z=(+Y)**2
Z
Out[14]: array([[ 1.
,
1.
,
[ 1.08329863,
1.08329863,
[ 1.1699292 ,
1.1699292 ,
...,
[ 8.51686797,
8.51686797,
[ 8.75676801,
8.75676801,
[ 9.
,
9.
,

1.
, 1.
, ...,
1.
],
1.08329863, 1.08329863, ...,
1.08329863],
1.1699292 , 1.1699292 , ...,
1.1699292 ],
8.51686797, 8.51686797, ...,
8.51686797],
8.75676801, 8.75676801, ...,
8.75676801],
9.
, 9.
, ...,
9.
]])
2

1.

1.08329863,
1.1699292 ,

8.51686797,
8.75676801,
9.

In [15]: Z[0:2]
Out[15]: array([[ 1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
[ 1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,

1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,

1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,

1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,

1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
,
1.
],
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863,
1.08329863]])

In [16]: Z[0:-1]
Out[16]: array([[ 1.
,
1.
,
[ 1.08329863,
1.08329863,
[ 1.1699292 ,
1.1699292 ,
...,
[ 8.28029988,
8.28029988,
[ 8.51686797,
8.51686797,
[ 8.75676801,
8.75676801,

1.
, 1.
, ...,
1.
],
1.08329863, 1.08329863, ...,
1.08329863],
1.1699292 , 1.1699292 , ...,
1.1699292 ],

1.

1.08329863,
1.1699292 ,

8.28029988, 8.28029988, ...,

8.28029988],
8.51686797, 8.51686797, ...,
8.51686797],
8.75676801, 8.75676801, ...,
8.75676801]])

8.28029988,

1.25989171, 1.25989171, ...,

1.25989171],
1.35318617, 1.35318617, ...,
1.35318617],
1.44981258, 1.44981258, ...,
1.44981258],

1.25989171,

8.51686797,
8.75676801,

In [17]: Z[3:]
Out[17]: array([[ 1.25989171,
1.25989171,
[ 1.35318617,
1.35318617,
[ 1.44981258,
1.44981258,
...,
[ 8.51686797,

8.51686797,
3

8.51686797, ...,

1.35318617,
1.44981258,

8.51686797,

8.51686797,
[ 8.75676801,
8.75676801,
[ 9.
,
9.
,

8.51686797],
8.75676801, 8.75676801, ...,
8.75676801],
9.
, 9.
, ...,
9.
]])

8.75676801,
9.

In [18]: Z[::-50]
Out[18]: array([[ 9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

In [19]: f=lambda a,b:a+10*b

F=fromfunction(f,[6,6])
F
Out[19]: array([[
[
[
[
[
[

0.,
1.,
2.,
3.,
4.,
5.,

10.,
11.,
12.,
13.,
14.,
15.,

20.,
21.,
22.,
23.,
24.,
25.,

30.,
31.,
32.,
33.,
34.,
35.,

40.,
41.,
42.,
43.,
44.,
45.,

4.,

5.])

50.],
51.],
52.],
53.],
54.],
55.]])

In [20]: F[:,0]
Out[20]: array([ 0.,

1.,

2.,

3.,

In [21]: F[:,3]
Out[21]: array([ 30.,

31.,

32.,

33.,

34.,

35.])

10.,

20.,

30.,

40.,

50.])

In [22]: F[0,:]
Out[22]: array([

0.,

In [23]: F[:3,:3]
Out[23]: array([[
[
[

0.,
1.,
2.,

10.,
11.,
12.,

20.],
21.],
22.]])

In [24]: F[:,0]
Out[24]: array([ 0.,

1.,

2.,

3.,

4.,

In [65]: F[1:4,2:4]
Out[65]: array([[ 21.,
[ 22.,
[ 23.,

31.],
32.],
33.]])
4

5.])

9.,
9.,
9.,
9.,

9.,
9.,
9.,
9.,

9., 9.,
9., 9.,
9., 9.,
9.]])

9.,
9.,
9.,

We can access data by issueing n:m which means from element in row and then call ,n:m to limit
it in specific column
In [26]: F[::-1]
Out[26]: array([[
[
[
[
[
[

5.,
4.,
3.,
2.,
1.,
0.,

15.,
14.,
13.,
12.,
11.,
10.,

25.,
24.,
23.,
22.,
21.,
20.,

35.,
34.,
33.,
32.,
31.,
30.,

45.,
44.,
43.,
42.,
41.,
40.,

55.],
54.],
53.],
52.],
51.],
50.]])

5.,
3.,
1.,

15.,
13.,
11.,

25.,
23.,
21.,

35.,
33.,
31.,

45.,
43.,
41.,

55.],
53.],
51.]])

In [27]: F[::-2]
Out[27]: array([[
[
[

In [28]: F[3::2,1::3]
Out[28]: array([[ 13.,
[ 15.,

43.],
45.]])

In [29]: r,h=symbols('r h')

In [30]: Area=2*pi*r**2 +2*pi*r*h
Area
Out[30]:
2hr + 2r2
In [31]: Volume=pi*r**2 *h
Volume
Out[31]:
hr2
In [32]: h_r=solve(Volume-1)[0]
h_r
Out[32]:


1
h: 2
r

In [33]: Area_r=Area.subs(h_r)
Area_r

Out[33]:
2r2 +

2
r

In [34]: rsol=solve(Area_r.diff(r))[0]
rsol
Out[34]:
2

23

23
In [35]: _.evalf()
Out[35]:
0.541926070139289
In [36]: Area_r.diff(r,2).subs(r,rsol)
Out[36]:
12
In [37]: Area_r.subs(r,rsol)
Out[37]:

3
3 23
In [38]: _.evalf()
Out[38]:
5.53581044593209
In [39]: func_r=lambdify(r,2*pi*r**2 +2/r)
func_r(r)
Out[39]:
6.28318530717959r2 +
In [68]: plotting.plot(func_r(r),(r,0,2))
<IPython.core.display.Javascript object>

<IPython.core.display.HTML object>

2
r

Out[68]: <sympy.plotting.plot.Plot at 0x7f76daf77898>

In [41]: from scipy.optimize import *
In [67]: min_r=minimize_scalar(func_r,bracket=[0.1,4])
min_r
Out[67]:

0.0.1

fun:
nfev:
nit:
success:
x:

5.5358104459320856
16
15
True
0.54192607725571351

Newtons method

In [43]: x=symbols('x')
x_k1,x_k,k=symbols('x_k1 x_k k')
f=Function('f')
x_k1=x_k -f(k).diff(k)/f(k).diff(k,2)
x_k1
Out[43]:
xk

0.0.2

d
dk f (k)
d2
f (k)
dk2

Newtons method for multivariate question

It should looks like this!

x_k1=x_k - hessian(f,x_k)**-1 *gradient(f)
In [44]: x1,x2=symbols('x_1 x_2')
In [45]: f_sym=(x1-1)**4 +5*(x2-1)**2 -2*x1*x2
f_sym
Out[45]:
2x1 x2 + (x1 1)4 + 5 (x2 1)2
In [46]: #Gradient
fprime_sym=Matrix([f_sym]).jacobian([x1,x2]).tolist()[0]
#fprime_sym=[f_sym.diff(x_) for x_ in (x1, x2)]
fprime_sym
Out[46]:
h
2x2 + 4 (x1 1)3 ,

i
2x1 + 10x2 10

In [47]: #Hessian
fhess_sym=hessian(f_sym,[x1,x2]).tolist()
#fhess_sym = [[f_sym.diff(x1_, x2_) for x1_ in (x1, x2)] for x2_ in (x1, x
fhess_sym
Out[47]:
hh

12 (x1 1)2 ,

i
2 ,

[2,

i
10]

In [48]: f_lmbda=lambdify([x1,x2],f_sym)
f_lmbda(x1,x2)
Out[48]:
2x1 x2 + (x1 1)4 + 5 (x2 1)2
In [49]: fprime_lmbda=lambdify([x1,x2],fprime_sym)
fprime_lmbda(x1,x2)
Out[49]:
h

2x2 + 4 (x1 1)3 ,

i
2x1 + 10x2 10

In [50]: fhess_lambda=lambdify([x1,x2],fhess_sym)
fhess_lambda(x1,x2)
Out[50]:
hh

12 (x1 1)2 ,

i
2 ,

[2,

i
10]

In [51]: def func_XY_to_X_Y(f):

"""
Wrapper for f(X) -> f(X[0], X[1])
"""
return lambda X:array(f(X[0],X[1]))
In [52]: f=func_XY_to_X_Y(f_lmbda)
fprime=func_XY_to_X_Y(fprime_lmbda)
fhess=func_XY_to_X_Y(fhess_lambda)
In [53]: x_opt=fmin_ncg(f,[0,0],fprime=fprime,fhess=fhess)
x_opt1=optimize.fmin_bfgs(f, (0, 0))
Optimization terminated successfully.
Current function value: -3.867223
Iterations: 8
Function evaluations: 10
Gradient evaluations: 17
Hessian evaluations: 8
8

Optimization terminated successfully.

Current function value: -3.867223
Iterations: 10
Function evaluations: 56
Gradient evaluations: 14

In [54]: x_opt.tolist(),x_opt1.tolist()
Out[54]:

([1.882926129296321,

1.3765852258592641] ,

[1.8829260403989003,

1.3765852205473148])

In [69]: fig, ax = plt.subplots(figsize=(6, 4))

x_ = y_ = linspace(-1, 4, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, f_lmbda(X, Y), 50)
ax.plot(x_opt[0], x_opt[1], 'r*', markersize=15)
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
plt.colorbar(c, ax=ax)
show()
<IPython.core.display.Javascript object>

<IPython.core.display.HTML object>

In [56]: fig, ax = plt.subplots(figsize=(6, 4))

x_ = y_ = linspace(-1, 4, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, f_lmbda(X, Y), 50)
ax.plot(x_opt1[0], x_opt1[1], 'r*', markersize=15)
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
plt.colorbar(c, ax=ax)
show()
<IPython.core.display.Javascript object>

<IPython.core.display.HTML object>

In [57]: x=x0,x1,x2,l=symbols('x_0 x_1 x_2 lambda')

x=[x0,x1,x2,l]
f=x0*x1*x2
f
9

Out[57]:
x0 x1 x2
In [58]: g=2*(x0*x1+x1*x2+x0*x2)-1
g
Out[58]:
2x0 x1 + 2x0 x2 + 2x1 x2 1
In [59]: L=f+l*g
L
Out[59]:
(2x0 x1 + 2x0 x2 + 2x1 x2 1) + x0 x1 x2
In [60]: grad_L=Matrix([L]).jacobian(x)
#grad_L = [diff(L, x_) for x_ in x]
grad_L
Out[60]:


(2x1 + 2x2 ) + x1 x2 (2x0 + 2x2 ) + x0 x2 (2x0 + 2x1 ) + x0 x1 2x0 x1 + 2x0 x2 + 2x1 x2 1
In [61]: sols=solve(grad_L)
sols
Out[61]:
"(

6
,
:
24

6
x0 :
,
6

6
x1 :
,
6

)
6
x2 :
,
6

6
:
,
24

6
x0 :
,
6

We can use *X to unpack the value in the given array or list.

In [62]: x=1,2,3,4,5
In [ ]:

10

6
x1 :
,
6

)#
6
x2 :
6