Changeset 902

Show
Ignore:
Timestamp:
09/23/05 11:14:37 (3 years ago)
Author:
rkern
Message:

r3087@Blasphemy: kern | 2005-09-23 09:13:29 -0700
More tutorial stuff.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • nbdoc/trunk/scipy_tutorial/chapter01.pybk

    r901 r902  
    1 <notebook version="1"><head/><ipython-log id="default-log"/><sheet id="chapter01"> 
     1<notebook version="1"><head/><ipython-log id="default-log"><cell number="0"><input>   
     2############DO NOT EDIT THIS CELL############   
     3from pylab import *   
     4switch_backend('WXAgg')   
     5ion()   
     6</input></cell><cell number="1"><input> 
     7from scipy import * 
     8</input></cell><cell number="2"><input> 
     9info(optimize.fmin) 
     10</input></cell><cell number="3"><input> 
     11optimize.fmin? 
     12</input></cell></ipython-log><sheet id="chapter01"> 
    213<title>Introduction</title> 
    3 <para>SciPy is a collection of mathematical algorithms and convenience functions 
    4 built on the Numeric extension for Python. It adds significant power to the 
    5 interactive Python session by exposing the user to high-level commands and 
    6 classes for the manipulation and visualization of data. With SciPy, an 
    7 interactive Python session becomes a data-processing and system-prototyping 
    8 environment rivaling sytems such as Matlab, IDL, Octave, R-Lab, and SciLab.  
    9 </para> 
    10 <para>The additional power of using SciPy within Python, however, is that a powerful 
    11 programming language is also available for use in developing sophisticated 
    12 programs and specialized applications. Scientific applications written in SciPy 
    13 benefit from the development of additional modules in numerous niche's of the 
    14 software landscape by developers across the world. Everything from parallel 
    15 programming to web and data-base subroutines and classes have been made 
    16 available to the Python programmer. All of this power is available in addition 
    17 to the mathematical libraries in SciPy. 
    18 </para> 
    19 <para>This document provides a tutorial for the first-time user of SciPy to help get 
    20 started with some of the features available in this powerful package. It is 
    21 assumed that the user has already installed the package. Some general Python 
    22 facility is also assumed such as could be acquired by working through the 
    23 Tutorial in the Python distribution. For further introductory help the user is 
    24 directed to the Numeric documentation. Throughout this tutorial it is assumed 
    25 that the user has imported all of the names defined in the SciPy namespace using 
    26 the command 
    27 </para> 
    28 <programlisting>XXX: ipython 
    29 >>> from scipy import * 
    30 </programlisting> 
    31  
     14<para>SciPy is a collection of mathematical algorithms and convenience functions built on the Numeric extension for Python. It adds significant power to the interactive Python session by exposing the user to high-level commands and classes for the manipulation and visualization of data. With SciPy, an interactive Python session becomes a data-processing and system-prototyping environment rivaling sytems such as Matlab, IDL, Octave, R-Lab, and SciLab. </para><para>The additional power of using SciPy within Python, however, is that a powerful programming language is also available for use in developing sophisticated programs and specialized applications. Scientific applications written in SciPy benefit from the development of additional modules in numerous niche's of the software landscape by developers across the world. Everything from parallel programming to web and data-base subroutines and classes have been made available to the Python programmer. All of this power is available in addition to the mathematical libraries in SciPy. </para><para>This document provides a tutorial for the first-time user of SciPy to help get started with some of the features available in this powerful package. It is assumed that the user has already installed the package. Some general Python facility is also assumed such as could be acquired by working through the Tutorial in the Python distribution. For further introductory help the user is directed to the Numeric documentation. Throughout this tutorial it is assumed that the user has imported all of the names defined in the SciPy namespace using the command </para><ipython-block logid="default-log"><ipython-cell type="input" number="1"/></ipython-block> 
    3215<section> 
    3316<title>General Help</title> 
    …… 
    5134in that module is printed. For example: 
    5235</para> 
    53 <programlisting>XXX: ipython 
    54 >>> info(optimize.fmin) 
    55  fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, 
    56       full_output=0, printmessg=1) 
    57  
    58 Minimize a function using the simplex algorithm. 
    59  
    60 Description: 
    61  
    62   Uses a Nelder-Mead simplex algorithm to find the minimum of function 
    63   of one or more variables. 
    64  
    65 Inputs: 
    66  
    67   func -- the Python function or method to be minimized. 
    68   x0 -- the initial guess. 
    69   args -- extra arguments for func. 
    70   xtol -- relative tolerance 
    71  
    72 Outputs: (xopt, {fopt, warnflag}) 
    73  
    74   xopt -- minimizer of function 
    75  
    76   fopt -- value of function at minimum: fopt = func(xopt) 
    77   warnflag -- Integer warning flag: 
    78               1 : 'Maximum number of function evaluations.' 
    79               2 : 'Maximum number of iterations.' 
    80  
    81 Additional Inputs: 
    82  
    83   xtol -- acceptable relative error in xopt for convergence. 
    84   ftol -- acceptable relative error in func(xopt) for convergence. 
    85   maxiter -- the maximum number of iterations to perform. 
    86   maxfun -- the maximum number of function evaluations. 
    87   full_output -- non-zero if fval and warnflag outputs are desired. 
    88   printmessg -- non-zero to print convergence messages. 
    89    
    90 </programlisting> 
     36<ipython-block logid="default-log"><ipython-cell type="input" number="2"/></ipython-block><para>(For some reason nbshell does not capture the output of this function, yet. The implementation of info() is very naughty and uses sys.stdout as a default argument instead of using None, then setting the default value of sys.stdout inside the function. If that had been the case, then we would have captured the output. </para><para>Naughty scipy!) </para> 
    9137<para>Another useful command is <code>source</code>. When given a function written in 
    9238Python as an argument, it prints out a listing of the source code for that 
    …… 
    10753<code>scipy_base</code> package is provided as well.  
    10854</para> 
    109 <para>Two other packages are installed at the higher-level: 
    110 scipy_distutils and weave. These two packages while distributed with main scipy 
    111 package could see use independently of scipy and so are treated as separate 
    112 packages and described elsewhere. 
    113 </para> 
    114 <para>The remaining subpackages are summarized in the following table (a * denotes an 
    115 optional sub-package that requires additional libraries to function or is not 
    116 available on all platforms).  
    117 </para> 
    118 <para>XXX: Table here.</para> 
    119 <para>Because of their ubiquitousness, some of the functions in these 
    120 subpackages are also made available in the scipy namespace to ease their use in 
    121 interactive sessions and programs. In addition, many convenience functions are 
    122 located in the scipy_base package and the in the top-level of the scipy package. 
    123 Before looking at the sub-packages individually, we will first look at some of 
    124 these common functions.  
    125 </para> 
    126 </section> 
     55<para>Two other packages are installed at the higher-level: scipy_distutils and weave. These two packages while distributed with main scipy package could see use independently of scipy and so are treated as separate packages and described elsewhere. </para><para>The remaining subpackages are summarized in the following table (a * denotes an optional sub-package that requires additional libraries to function or is not available on all platforms). </para><para>XXX: Table here. </para><para>Because of their ubiquitousness, some of the functions in these subpackages are also made available in the scipy namespace to ease their use in interactive sessions and programs. In addition, many convenience functions are located in the scipy_base package and the in the top-level of the scipy package. Before looking at the sub-packages individually, we will first look at some of these common functions. </para></section> 
    12756</sheet></notebook> 

  • nbdoc/trunk/scipy_tutorial/chapter02.pybk

    r901 r902  
    1 <notebook version="1"><head/><ipython-log id="default-log"/><sheet id="chapter02"> 
     1<notebook version="1"><head/><ipython-log id="default-log"><cell number="0"><input>   
     2############DO NOT EDIT THIS CELL############   
     3from pylab import *   
     4switch_backend('WXAgg')   
     5ion()   
     6</input></cell><cell number="1"><input> 
     7import scipy; from scipy import * 
     8</input></cell><cell number="2"><input> 
     9scipy.sqrt(-1) 
     10</input><output> 
     111j 
     12</output></cell><cell number="3"><input> 
     13concatenate(([3], [0]*5, arange(-1,1.002,2/9.0)) 
     14) 
     15</input><output> 
     16[ 3.        , 0.        , 0.        , 0.        , 0.        , 0. 
     17, 
     18      -1.        ,-0.77777778,-0.55555556,-0.33333333,-0.11111111, 
     190.11111111, 
     20       0.33333333, 0.55555556, 0.77777778, 1.        ,] 
     21</output></cell><cell number="4"><input> 
     22r_[3, [0]*5, -1:1:10j] 
     23</input><output> 
     24[ 3.        , 0.        , 0.        , 0.        , 0.        , 0. 
     25, 
     26      -1.        ,-0.77777778,-0.55555556,-0.33333333,-0.11111111, 
     270.11111111, 
     28       0.33333333, 0.55555556, 0.77777778, 1.        ,] 
     29</output></cell><cell number="5"><input> 
     30mgrid[0:5, 0:5] 
     31</input><output> 
     32[[[0,0,0,0,0,] 
     33  [1,1,1,1,1,] 
     34  [2,2,2,2,2,] 
     35  [3,3,3,3,3,] 
     36  [4,4,4,4,4,]] 
     37 [[0,1,2,3,4,] 
     38  [0,1,2,3,4,] 
     39  [0,1,2,3,4,] 
     40  [0,1,2,3,4,] 
     41  [0,1,2,3,4,]]] 
     42</output></cell><cell number="6"><input> 
     43mgrid[0:5:4j, 0:5:4j] 
     44</input><output> 
     45[[[ 0.        , 0.        , 0.        , 0.        ,] 
     46  [ 1.66666667, 1.66666667, 1.66666667, 1.66666667,] 
     47  [ 3.33333333, 3.33333333, 3.33333333, 3.33333333,] 
     48  [ 5.        , 5.        , 5.        , 5.        ,]] 
     49 [[ 0.        , 1.66666667, 3.33333333, 5.        ,] 
     50  [ 0.        , 1.66666667, 3.33333333, 5.        ,] 
     51  [ 0.        , 1.66666667, 3.33333333, 5.        ,] 
     52  [ 0.        , 1.66666667, 3.33333333, 5.        ,]]] 
     53</output></cell><cell number="7"><input> 
     54p = poly1d([3,4,5]) 
     55</input></cell><cell number="8"><input> 
     56print p 
     57</input><stdout> 
     58   2 
     593 x + 4 x + 5 
     60</stdout></cell><cell number="9"><input> 
     61print p*p 
     62</input><stdout> 
     63   4      3      2 
     649 x + 24 x + 46 x + 40 x + 25 
     65</stdout></cell><cell number="10"><input> 
     66print p.integ(k=6) 
     67</input><stdout> 
     68 3     2 
     69x + 2 x + 5 x + 6 
     70</stdout></cell><cell number="11"><input> 
     71print p.deriv() 
     72</input><stdout> 
     73 
     746 x + 4 
     75</stdout></cell><cell number="12"><input> 
     76p([4,5]) 
     77</input><output> 
     78[ 69,100,] 
     79</output></cell><cell number="13"><input> 
     80def addsubtract(a, b): 
     81    if a &gt; b: 
     82        return a - b 
     83    else: 
     84        return a + b 
     85         
     86</input></cell><cell number="14"><input> 
     87vec_addsubtract = vectorize(addsubtract) 
     88</input></cell><cell number="15"><input> 
     89vec_addsubtract([0,3,6,9], [1,3,5,7]) 
     90</input><output> 
     91[1,6,1,2,] 
     92</output></cell><cell number="16"><input> 
     93x = r_[-2:3] 
     94</input></cell><cell number="17"><input> 
     95x 
     96</input><output> 
     97[-2,-1, 0, 1, 2,] 
     98</output></cell><cell number="19"><input> 
     99select([x &gt; 3, x &gt;= 0], [0,x+2]) 
     100</input><output> 
     101[0,0,2,3,4,] 
     102</output></cell></ipython-log><sheet id="chapter02"> 
    2103<title>Basic functions in scipy_base and top-level scipy</title> 
    3104<section> 
    …… 
    7108additionally importing Numeric. In addition, the universal functions (addition, 
    8109subtraction, division) have been altered to not raise exceptions if 
    9 floating-point errors are encountered<footnote><para>These changes are all made in a new module 
    10 (fastumath) that is part of the scipy_base package. The old functionality is 
    11 still available in umath (part of Numeric) if you need it (note: importing umath 
    12 or fastumath resets the behavior of the infix operators to use the umath or 
    13 fastumath ufuncs respectively).</para></footnote>, instead NaN's and Inf's are returned in the 
     110floating-point errors are encountered<footnote><para>These changes are all made in a new module (fastumath) that is part of the scipy_base package. The old functionality is still available in umath (part of Numeric) if you need it (note: importing umath or fastumath resets the behavior of the infix operators to use the umath or fastumath ufuncs respectively). </para></footnote>, instead NaN's and Inf's are returned in the 
    14111arrays. To assist in detection of these events new universal functions (isnan, 
    15112isfinite, isinf) have been added. 
    …… 
    20117operations (except logical_XXX functions) all return arrays of unsigned bytes 
    21118(8-bits per element instead of the old 32-bits, or even 64-bits) per 
    22 element<footnote><para>Be 
    23 careful when treating logical expressions as integers as the 8-bit integers may 
    24 silently overflow at 256.</para></footnote>. 
     119element<footnote><para>Be careful when treating logical expressions as integers as the 8-bit integers may silently overflow at 256. </para></footnote>. 
    25120</para> 
    26 <para>In an effort to get a consistency for default arguments, some of the 
    27 default arguments have changed from Numeric. The idea is for you to use scipy as 
    28 a base package instead of Numeric anyway. 
    29 </para> 
    30 <para>Finally, some of the basic functions like log, sqrt, and inverse trig functions 
    31 have been modified to return complex numbers instead of NaN's where appropriate. 
    32 </para> 
    33 <programlisting>XXX: ipython 
    34 >>> scipy.sqrt(-1) </programlisting> 
     121<para>In an effort to get a consistency for default arguments, some of the default arguments have changed from Numeric. The idea is for you to use scipy as a base package instead of Numeric anyway. </para><para>Finally, some of the basic functions like log, sqrt, and inverse trig functions have been modified to return complex numbers instead of NaN's where appropriate. </para><ipython-block logid="default-log"><ipython-cell type="input" number="1"/><ipython-cell type="input" number="2"/><ipython-cell type="output" number="2"/></ipython-block> 
    35122</section> 
    36123<section> 
    37124<title>Alter numeric</title> 
    38 <para>With the command scipy.alter_numeric() you can now use index and mask arrays inside brackets and the coercion rules of Numeric will be changed so that Python scalars will not be used to determine the type of the output of an expression. 
    39 </para> 
    40 </section> 
     125<para>With the command scipy.alter_numeric() you can now use index and mask arrays inside brackets and the coercion rules of Numeric will be changed so that Python scalars will not be used to determine the type of the output of an expression. </para></section> 
    41126<section> 
    42127<title>scipy_base routines</title> 
    43 <para>The purpose of scipy_base is to collect general-purpose routines 
    44 that the other sub-packages can use and to provide a simple replacement for 
    45 Numeric. Anytime you might think to import Numeric, you can import scipy_base 
    46 instead and remove yourself from direct dependence on Numeric. These routines 
    47 are divided into several files for organizational purposes, but they are all 
    48 available under the scipy_base namespace (and the scipy namespace). There are 
    49 routines for type handling and type checking, shape and matrix manipulation, 
    50 polynomial processing, and other useful functions. Rather than giving a detailed 
    51 description of each of these functions (which is available using the help, info 
    52 and source commands), this tutorial will discuss some of the more useful 
    53 commands which require a little introduction to use to their full potential.  
    54 </para> 
    55 <section> 
     128<para>The purpose of scipy_base is to collect general-purpose routines that the other sub-packages can use and to provide a simple replacement for Numeric. Anytime you might think to import Numeric, you can import scipy_base instead and remove yourself from direct dependence on Numeric. These routines are divided into several files for organizational purposes, but they are all available under the scipy_base namespace (and the scipy namespace). There are routines for type handling and type checking, shape and matrix manipulation, polynomial processing, and other useful functions. Rather than giving a detailed description of each of these functions (which is available using the help, info and source commands), this tutorial will discuss some of the more useful commands which require a little introduction to use to their full potential. </para><section> 
    56129<title>Type handling</title> 
    57 <para>Note the difference between iscomplex (isreal) and iscomplexobj 
    58 (isrealobj). The former command is array based and returns byte arrays of ones 
    59 and zeros providing the result of the element-wise test. The latter command is 
    60 object based and returns a scalar describing the result of the test on the 
    61 entire object.  
    62 </para> 
    63 <para>Often it is required to get just the real and/or imaginary part of a complex 
    64 number. While complex numbers and arrays have attributes that return those 
    65 values, if one is not sure whether or not the object will be complex-valued, it 
    66 is better to use the functional forms real and imag. These functions succeed for 
    67 anything that can be turned into a Numeric array. Consider also the function 
    68 real_if_close which transforms a complex-valued number with tiny imaginary part 
    69 into a real number.  
    70 </para> 
    71 <para>Occasionally the need to check whether or not a number is a scalar (Python 
    72 (long)int, Python float, Python complex, or rank-0 array) occurs in coding. This 
    73 functionality is provided in the convenient function isscalar which returns a 1 
    74 or a 0.  
    75 </para> 
    76 <para>Finally, ensuring that objects are a certain Numeric type occurs often enough 
    77 that it has been given a convenient interface in SciPy through the use of the 
    78 cast dictionary. The dictionary is keyed by the type it is desired to cast to 
    79 and the dictionary stores functions to perform the casting. Thus, &gt;&gt;&gt; a = 
    80 cast['f'](d) returns an array of float32 from d. This function is also useful as 
    81 an easy way to get a scalar of a certain type: &gt;&gt;&gt; fpi = cast['f'](pi) although 
    82 this should not be needed if you use alter_numeric(). 
    83 </para> 
    84 </section> 
     130<para>Note the difference between iscomplex (isreal) and iscomplexobj (isrealobj). The former command is array based and returns byte arrays of ones and zeros providing the result of the element-wise test. The latter command is object based and returns a scalar describing the result of the test on the entire object. </para><para>Often it is required to get just the real and/or imaginary part of a complex number. While complex numbers and arrays have attributes that return those values, if one is not sure whether or not the object will be complex-valued, it is better to use the functional forms real and imag. These functions succeed for anything that can be turned into a Numeric array. Consider also the function real_if_close which transforms a complex-valued number with tiny imaginary part into a real number. </para><para>Occasionally the need to check whether or not a number is a scalar (Python (long)int, Python float, Python complex, or rank-0 array) occurs in coding. This functionality is provided in the convenient function isscalar which returns a 1 or a 0. </para><para>Finally, ensuring that objects are a certain Numeric type occurs often enough that it has been given a convenient interface in SciPy through the use of the cast dictionary. The dictionary is keyed by the type it is desired to cast to and the dictionary stores functions to perform the casting. Thus, &gt;&gt;&gt; a = cast['f'](d) returns an array of float32 from d. This function is also useful as an easy way to get a scalar of a certain type: &gt;&gt;&gt; fpi = cast['f'](pi) although this should not be needed if you use alter_numeric(). </para></section> 
    85131<section> 
    86132<title>Index tricks</title> 
    87 <para>There are some class instances that make special use of the slicing 
    88 functionality to provide efficient means for array construction. This part will 
    89 discuss the operation of mgrid, ogrid, r_, and c_ for quickly constructing 
    90 arrays.  
    91 </para> 
    92 <para>One familiar with Matlab may complain that it is difficult to construct arrays 
    93 from the interactive session with Python. Suppose, for example that one wants to 
    94 construct an array that begins with 3 followed by 5 zeros and then contains 10 
    95 numbers spanning the range -1 to 1 (inclusive on both ends). Before SciPy, you 
    96 would need to enter something like the following &gt;&gt;&gt; 
    97 concatenate(([3],[0]*5,arange(-1,1.002,2/9.0)). With the r_ command one can 
    98 enter this as &gt;&gt;&gt; r_[3,[0]*5,-1:1:10j] which can ease typing and make for more 
    99 readable code. Notice how objects are concatenated, and the slicing syntax is 
    100 (ab)used to construct ranges. The other term that deserves a little explanation 
    101 is the use of the complex number 10j as the step size in the slicing syntax. 
    102 This non-standard use allows the number to be interpreted as the number of 
    103 points to produce in the range rather than as a step size (note we would have 
    104 used the long integer notation, 10L, but this notation may go away in Python as 
    105 the integers become unified). This non-standard usage may be unsightly to some, 
    106 but it gives the user the ability to quickly construct complicated vectors in a 
    107 very readable fashion. When the number of points is specified in this way, the 
    108 end-point is inclusive. 
    109 </para> 
    110 <para>The "r" stands for row concatenation because if the objects between commas are 2 
    111 dimensional arrays, they are stacked by rows (and thus must have commensurate 
    112 columns). There is an equivalent command c_ that stacks 2d arrays by columns but 
    113 works identically to r_ for 1d arrays.  
    114 </para> 
    115 <para>Another very useful class instance which makes use of extended slicing notation 
    116 is the function mgrid. In the simplest case, this function can be used to 
    117 construct 1d ranges as a convenient substitute for arange. It also allows the 
    118 use of complex-numbers in the step-size to indicate the number of points to 
    119 place between the (inclusive) end-points. The real purpose of this function 
    120 however is to produce N, N-d arrays which provide coordinate arrays for an 
    121 N-dimensional volume. The easiest way to understand this is with an example of 
    122 its usage:  
    123 </para> 
    124 <programlisting>XXX:  
    125 &gt;&gt;&gt; mgrid[0:5,0:5] 
    126 array([[[0, 0, 0, 0, 0], 
    127         [1, 1, 1, 1, 1], 
    128         [2, 2, 2, 2, 2], 
    129         [3, 3, 3, 3, 3], 
    130         [4, 4, 4, 4, 4]], 
    131        [[0, 1, 2, 3, 4], 
    132         [0, 1, 2, 3, 4], 
    133         [0, 1, 2, 3, 4], 
    134         [0, 1, 2, 3, 4], 
    135         [0, 1, 2, 3, 4]]])   
    136 &gt;&gt;&gt; mgrid[0:5:4j,0:5:4j] 
    137 array([[[ 0.    ,  0.    ,  0.    ,  0.    ], 
    138         [ 1.6667,  1.6667,  1.6667,  1.6667], 
    139         [ 3.3333,  3.3333,  3.3333,  3.3333], 
    140         [ 5.    ,  5.    ,  5.    ,  5.    ]], 
    141        [[ 0.    ,  1.6667,  3.3333,  5.    ], 
    142         [ 0.    ,  1.6667,  3.3333,  5.    ], 
    143         [ 0.    ,  1.6667,  3.3333,  5.    ], 
    144         [ 0.    ,  1.6667,  3.3333,  5.    ]]]) 
    145 </programlisting> 
    146 <para>Having meshed arrays like this is sometimes very useful. However, it is 
    147 not always needed just to evaluate some N-dimensional function over a grid due 
    148 to the array-broadcasting rules of Numeric and SciPy. If this is the only 
    149 purpose for generating a meshgrid, you should instead use the function ogrid 
    150 which generates an "open" grid using NewAxis judiciously to create N, N-d arrays 
    151 where only one-dimension in each array has length greater than 1. This will save 
    152 memory and create the same result if the only purpose for the meshgrid is to 
    153 generate sample points for evaluation of an N-d function.  
    154 </para> 
    155 </section> 
     133<para>There are some class instances that make special use of the slicing functionality to provide efficient means for array construction. This part will discuss the operation of mgrid, ogrid, r_, and c_ for quickly constructing arrays. </para><para>One familiar with Matlab may complain that it is difficult to construct arrays from the interactive session with Python. Suppose, for example that one wants to construct an array that begins with 3 followed by 5 zeros and then contains 10 numbers spanning the range -1 to 1 (inclusive on both ends). Before SciPy, you would need to enter something like the following </para><ipython-block logid="default-log"><ipython-cell type="input" number="3"/><ipython-cell type="output" number="3"/></ipython-block><para>With the r_ command one can enter this as </para><ipython-block logid="default-log"><ipython-cell type="input" number="4"/><ipython-cell type="output" number="4"/></ipython-block><para>which can ease typing and make for more readable code. Notice how objects are concatenated, and the slicing syntax is (ab)used to construct ranges. The other term that deserves a little explanation is the use of the complex number 10j as the step size in the slicing syntax. This non-standard use allows the number to be interpreted as the number of points to produce in the range rather than as a step size (note we would have used the long integer notation, 10L, but this notation may go away in Python as the integers become unified). This non-standard usage may be unsightly to some, but it gives the user the ability to quickly construct complicated vectors in a very readable fashion. When the number of points is specified in this way, the end-point is inclusive. </para><para>The "r" stands for row concatenation because if the objects between commas are 2 dimensional arrays, they are stacked by rows (and thus must have commensurate columns). There is an equivalent command c_ that stacks 2d arrays by columns but works identically to r_ for 1d arrays. </para><para>Another very useful class instance which makes use of extended slicing notation is the function mgrid. In the simplest case, this function can be used to construct 1d ranges as a convenient substitute for arange. It also allows the use of complex-numbers in the step-size to indicate the number of points to place between the (inclusive) end- points. The real purpose of this function however is to produce N, N-d arrays which provide coordinate arrays for an N-dimensional volume. The easiest way to understand this is with an example of its usage: </para><ipython-block logid="default-log"><ipython-cell type="input" number="5"/><ipython-cell type="output" number="5"/><ipython-cell type="input" number="6"/><ipython-cell type="output" number="6"/></ipython-block> 
     134<para>Having meshed arrays like this is sometimes very useful. However, it is not always needed just to evaluate some N-dimensional function over a grid due to the array-broadcasting rules of Numeric and SciPy. If this is the only purpose for generating a meshgrid, you should instead use the function ogrid which generates an "open" grid using NewAxis judiciously to create N, N-d arrays where only one-dimension in each array has length greater than 1. This will save memory and create the same result if the only purpose for the meshgrid is to generate sample points for evaluation of an N-d function. </para></section> 
    156135 
    157136<section> 
    158137<title>Shape manipulation</title> 
    159 <para>In this category of functions are routines for squeezing out length-one 
    160 dimensions from N-dimensional arrays, ensuring that an array is at least 1-, 2-, 
    161 or 3-dimensional, and stacking (concatenating) arrays by rows, columns, and 
    162 "pages" (in the third dimension). Routines for splitting arrays (roughly the 
    163 opposite of stacking arrays) are also available.  
    164 </para> 
    165 </section> 
     138<para>In this category of functions are routines for squeezing out length- one dimensions from N-dimensional arrays, ensuring that an array is at least 1-, 2-, or 3-dimensional, and stacking (concatenating) arrays by rows, columns, and "pages" (in the third dimension). Routines for splitting arrays (roughly the opposite of stacking arrays) are also available. </para></section> 
    166139<section> 
    167140<title>Matrix manipulations</title> 
    168 <para>These are functions specifically suited for 2-dimensional arrays that were 
    169 part of MLab in the Numeric distribution, but have been placed in scipy_base for 
    170 completeness so that users are not importing Numeric.  
    171 </para> 
    172 </section> 
     141<para>These are functions specifically suited for 2-dimensional arrays that were part of MLab in the Numeric distribution, but have been placed in scipy_base for completeness so that users are not importing Numeric. </para></section> 
    173142<section> 
    174143<title>Polynomials</title> 
    175 <para>There are two (interchangeable) ways to deal with 1-d polynomials in 
    176 SciPy. The first is to use the poly1d class in scipy_base. This class accepts 
    177 coefficients or polynomial roots to initialize a polynomial. The polynomial 
    178 object can then be manipulated in algebraic expressions, integrated, 
    179 differentiated, and evaluated. It even prints like a polynomial: 
    180 </para> 
    181 <programlisting> XXX: ipython 
    182 &gt;&gt;&gt; p = 
    183 poly1d([3,4,5]) 
    184 &gt;&gt;&gt; print p 
    185    2 
    186 3 x + 4 x + 5 
    187 &gt;&gt;&gt; print p*p 
    188    4      3      2 
    189 9 x + 24 x + 46 x + 40 x + 25 
    190 &gt;&gt;&gt; print p.integ(k=6) 
    191  3     2 
    192 x + 2 x + 5 x + 6 
    193 &gt;&gt;&gt; print p.deriv() 
    194 6 x + 4 
    195 &gt;&gt;&gt; p([4,5]) 
    196 array([ 69, 100]) 
    197 </programlisting> 
    198 <para>The other way to handle polynomials is as an array of coefficients with 
    199 the first element of the array giving the coefficient of the highest power. 
    200 There are explicit functions to add, subtract, multiply, divide, integrate, 
    201 differentiate, and evaluate polynomials represented as sequences of 
    202 coefficients.  
    203 </para> 
    204 </section> 
     144<para>There are two (interchangeable) ways to deal with 1-d polynomials in SciPy. The first is to use the poly1d class in scipy_base. This class accepts coefficients or polynomial roots to initialize a polynomial. The polynomial object can then be manipulated in algebraic expressions, integrated, differentiated, and evaluated. It even prints like a polynomial: </para><ipython-block logid="default-log"><ipython-cell type="input" number="7"/><ipython-cell type="input" number="8"/><ipython-cell type="stdout" number="8"/><ipython-cell type="input" number="9"/><ipython-cell type="stdout" number="9"/><ipython-cell type="input" number="10"/><ipython-cell type="stdout" number="10"/><ipython-cell type="input" number="11"/><ipython-cell type="stdout" number="11"/><ipython-cell type="input" number="12"/><ipython-cell type="output" number="12"/></ipython-block> 
     145<para>The other way to handle polynomials is as an array of coefficients with the first element of the array giving the coefficient of the highest power. There are explicit functions to add, subtract, multiply, divide, integrate, differentiate, and evaluate polynomials represented as sequences of coefficients. </para></section> 
    205146 
    206147<section> 
    207148<title>Vectorizing functions (vectorize)</title> 
    208 <para>One of the features that SciPy provides is a class vectorize to convert an 
    209 ordinary Python function which accepts scalars and returns scalars into a 
    210 "vectorized-function" with the same broadcasting rules as other Numeric 
    211 functions (i.e. the Universal functions, or ufuncs). For example, suppose you 
    212 have a Python function named addsubtract defined as: 
    213 </para> 
    214 <programlisting>XXX: ipython 
    215 &gt;&gt;&gt; def addsubtract(a,b): 
    216     if a &gt; b: 
    217         return a - b 
    218     else: 
    219         return a + b 
    220 </programlisting> 
    221 <para>which defines a function of two scalar variables and returns a scalar 
    222 result. The class vectorize can be used to "vectorize" this function so that 
    223 returns a function which takes array arguments and returns an array result: 
    224 </para> 
    225 <programlisting> XXX: ipython 
    226 &gt;&gt;&gt; vec_addsubtract([0,3,6,9],[1,3,5,7]) 
    227 array([1, 6, 1, 2]) 
    228 </programlisting> 
    229 <para>This particular function could have been written in vector form without 
    230 the use of vectorize. But, what if the function you have written is the result 
    231 of some optimization or integration routine. Such functions can likely only be 
    232 vectorized using vectorize. 
    233 </para> 
    234 </section> 
     149<para>One of the features that SciPy provides is a class vectorize to convert an ordinary Python function which accepts scalars and returns scalars into a "vectorized-function" with the same broadcasting rules as other Numeric functions (i.e. the Universal functions, or ufuncs). For example, suppose you have a Python function named addsubtract defined as: </para><ipython-block logid="default-log"><ipython-cell type="input" number="13"/></ipython-block><para>which defines a function of two scalar variables and returns a scalar result. The class vectorize can be used to "vectorize" this function so that returns a function which takes array arguments and returns an array result: </para><ipython-block logid="default-log"><ipython-cell type="input" number="14"/><ipython-cell type="input" number="15"/><ipython-cell type="output" number="15"/></ipython-block> 
     150<para>This particular function could have been written in vector form without the use of vectorize. But, what if the function you have written is the result of some optimization or integration routine. Such functions can likely only be vectorized using vectorize. </para></section> 
    235151<section> 
    236152<title>Other useful functions</title> 
    237 <para>There are several other functions in the scipy_base package including most 
    238 of the other functions that are also in MLab that comes with the Numeric 
    239 package. The reason for duplicating these functions is to allow SciPy to 
    240 potentially alter their original interface and make it easier for users to know 
    241 how to get access to functions &gt;&gt;&gt; from scipy import *. 
    242 </para> 
    243 <para>New functions which should be mentioned are mod(x,y) which can replace x%y 
    244 when it is desired that the result take the sign of y instead of x. Also 
    245 included is fix which always rounds to the nearest integer towards zero. For 
    246 doing phase processing, the functions angle, and unwrap are also useful. Also, 
    247 the linspace and logspace functions return equally spaced samples in a linear or 
    248 log scale. Finally, mention should be made of the new function select which 
    249 extends the functionality of where to include multiple conditions and multiple 
    250 choices. The calling convention is select(condlist,choicelist,default=0). Select 
    251 is a vectorized form of the multiple if-statement. It allows rapid construction 
    252 of a function which returns an array of results based on a list of conditions. 
    253 Each element of the return array is taken from the array in a choicelist 
    254 corresponding to the first condition in condlist that is true. For example  
    255 </para> 
    256 <programlisting> XXX: ipython 
    257 &gt;&gt;&gt; x = r_[-2:3] 
    258 &gt;&gt;&gt; x 
    259 array([-2, -1,  0,  1,  2]) 
    260 &gt;&gt;&gt; select([x &gt; 3, x &gt;= 0],[0,x+2]) 
    261 array([0, 0, 2, 3, 4]) 
    262 </programlisting> 
     153<para>There are several other functions in the scipy_base package including most of the other functions that are also in MLab that comes with the Numeric package. The reason for duplicating these functions is to allow SciPy to potentially alter their original interface and make it easier for users to know how to get access to functions &gt;&gt;&gt; from scipy import *. </para><para>New functions which should be mentioned are mod(x,y) which can replace x%y when it is desired that the result take the sign of y instead of x. Also included is fix which always rounds to the nearest integer towards zero. For doing phase processing, the functions angle, and unwrap are also useful. Also, the linspace and logspace functions return equally spaced samples in a linear or log scale. Finally, mention should be made of the new function select which extends the functionality of where to include multiple conditions and multiple choices. The calling convention is select(condlist,choicelist,default=0). Select is a vectorized form of the multiple if-statement. It allows rapid construction of a function which returns an array of results based on a list of conditions. Each element of the return array is taken from the array in a choicelist corresponding to the first condition in condlist that is true. For example </para><ipython-block logid="default-log"><ipython-cell type="input" number="16"/><ipython-cell type="input" number="17"/><ipython-cell type="output" number="17"/><ipython-cell type="input" number="19"/><ipython-cell type="output" number="19"/></ipython-block> 
    263154</section> 
    264155</section> 

  • nbdoc/trunk/scipy_tutorial/chapter04.pybk

    r901 r902  
    1 <notebook version="1"><head/><ipython-log id="default-log"/><sheet id="chapter04"> 
     1<notebook version="1"><head/><ipython-log id="default-log"><cell number="0"><input>   
     2############DO NOT EDIT THIS CELL############   
     3from pylab import *   
     4switch_backend('WXAgg')   
     5ion()   
     6</input></cell><cell number="1"><input> 
     7import scipy; from scipy import * 
     8</input></cell><cell number="2"><input> 
     9help(integrate) 
     10</input><stdout> 
     11Help on package scipy.integrate in scipy: 
     12 
     13NAME 
     14    scipy.integrate 
     15 
     16FILE 
     17    /Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4 
     18/site-packages/scipy_complete-0.3.3_309.4626-py2.4-macosx-10.4-ppc.egg 
     19/scipy/integrate/__init__.py 
     20 
     21DESCRIPTION 
     22    Integration routines 
     23    ==================== 
     24 
     25     Methods for Integrating Functions given function object. 
     26 
     27       quad          -- General purpose integration. 
     28       dblquad       -- General purpose double integration. 
     29       tplquad       -- General purpose triple integration. 
     30       fixed_quad    -- Integrate func(x) using Gaussian quadrature of 
     31order n. 
     32       quadrature    -- Integrate with given tolerance using Gaussian 
     33quadrature. 
     34       romberg       -- Integrate func using Romberg integration. 
     35 
     36     Methods for Integrating Functions given fixed samples. 
     37 
     38       trapz         -- Use trapezoidal rule to compute integral from 
     39samples. 
     40       cumtrapz      -- Use trapezoidal rule to cumulatively compute 
     41integral. 
     42       simps         -- Use Simpson's rule to compute integral from 
     43samples. 
     44       romb          -- Use Romberg Integration to compute integral 
     45from 
     46                        (2**k + 1) evenly-spaced samples. 
     47 
     48       See the special module's orthogonal polynomials (special) for 
     49Gaussian 
     50          quadrature roots and weights for other weighting factors and 
     51regions. 
     52 
     53     Interface to numerical integrators of ODE systems. 
     54 
     55       odeint        -- General integration of ordinary differential 
     56equations. 
     57       ode           -- Integrate ODE using vode routine. 
     58 
     59PACKAGE CONTENTS 
     60    _odepack 
     61    _quadpack 
     62    common_routines 
     63    info_integrate 
     64    ode 
     65    odepack 
     66    quadpack 
     67    quadrature 
     68    setup_integrate 
     69    vode 
     70 
     71DATA 
     72    Inf = inf 
     73    inf = inf 
     74 
     75</stdout></cell><cell number="3"><input> 
     76result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5) 
     77</input></cell><cell number="4"><input> 
     78print result 
     79</input><stdout> 
     80(1.1178179380783251, 7.8663172163807182e-09) 
     81</stdout></cell><cell number="5"><input> 
     82I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5)-4.0/27*sqrt(2)*sin(4.5)+sqrt(2*pi)* 
     83special.fresnel(3/sqrt(pi))[0]) 
     84</input></cell><cell number="6"><input> 
     85print I 
     86 
     87</input><stdout> 
     881.11781793809 
     89</stdout></cell><cell number="7"><input> 
     90print abs(result[0] - I) 
     91</input><stdout> 
     921.03759223435e-11 
     93</stdout></cell><cell number="8"><input> 
     94from scipy.integrate import quad, Inf 
     95</input></cell><cell number="9"><input> 
     96def integrand(t, n, x): 
     97    return exp(-x*t)/t**n 
     98</input></cell><cell number="10"><input> 
     99def expint(n, x): 
     100    return quad(integrand, 1, Inf, args=(n,x))[0] 
     101</input></cell><cell number="11"><input> 
     102vec_expint = vectorize(expint) 
     103</input></cell><cell number="12"><input> 
     104vec_expint(3, arange(1.0, 4.0, 0.5)) 
     105</input><output> 
     106[ 0.10969197, 0.05673949, 0.03013338, 0.01629537, 0.00893065, 
     1070.00494538,] 
     108</output></cell><cell number="13"><input> 
     109special.expn(3, arange(1.0, 4.0, 0.5)) 
     110</input><output> 
     111[ 0.10969197, 0.05673949, 0.03013338, 0.01629537, 0.00893065, 
     1120.00494538,] 
     113</output></cell><cell number="14"><input> 
     114result = quad(lambda x: expint(3, x), 0, Inf) 
     115</input></cell><cell number="15"><input> 
     116print result 
     117</input><stdout> 
     118(0.33333333325010889, 2.8604070802797126e-09) 
     119</stdout></cell><cell number="16"><input> 
     120I3 = 1.0/3.0 
     121</input></cell><cell number="17"><input> 
     122print I3 
     123</input><stdout> 
     1240.333333333333 
     125</stdout></cell><cell number="18"><input> 
     126print I3 - result[0] 
     127</input><stdout> 
     1288.32244273496e-11 
     129</stdout></cell><cell number="19"><input> 
     130from scipy.integrate import dblquad 
     131</input></cell><cell number="20"><input> 
     132def I(n): 
     133  return dblquad(lambda t,x: exp(-x*t)/t**n, 0, Inf, lambda x: 1, lambda x: Inf) 
     134</input></cell><cell number="21"><input> 
     135print I(4) 
     136</input><stdout> 
     137(0.25000000000435768, 1.0518246209542075e-09) 
     138</stdout></cell><cell number="22"><input> 
     139print I(3) 
     140</input><stdout> 
     141(0.33333333325010889, 2.8604070802797126e-09) 
     142</stdout></cell><cell number="23"><input> 
     143print I(2) 
     144</input><stdout> 
     145(0.49999999999857531, 1.8855523253868967e-09) 
     146</stdout></cell><cell number="24"><input> 
     147from scipy.integrate import odeint 
     148</input></cell><cell number="25"><input> 
     149from scipy.special import gamma, airy 
     150</input></cell><cell number="26"><input> 
     151y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0) 
     152</input></cell><cell number="27"><input> 
     153y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0) 
     154</input></cell><cell number="28"><input> 
     155y0 = [y0_0, y1_0] 
     156</input></cell><cell number="29"><input> 
     157def func(y, t): 
     158  return [t*y[1], y[0]] 
     159</input></cell><cell number="30"><input> 
     160def gradient(y, t): 
     161  return [[0,t], [1,0]] 
     162   
     163   
     164</input></cell><cell number="31"><input> 
     165x = arange(0, 4.0, 0.01) 
     166</input></cell><cell number="32"><input> 
     167t = x 
     168</input></cell><cell number="33"><input> 
     169ychk = airy(x)[0] 
     170</input></cell><cell number="34"><input> 
     171y = odeint(func, y0, t) 
     172</input></cell><cell number="35"><input> 
     173y2 = odeint(func, y0, t, Dfun=gradient) 
     174</input></cell><cell number="36"><input> 
     175import sys 
     176</input></cell><cell number="37"><input> 
     177sys.float_output_precision = 6 
     178</input></cell><cell number="38"><input> 
     179print ychk[:36:6] 
     180</input><stdout> 
     181[ 0.355028, 0.339511, 0.324068, 0.308763, 0.293658, 0.278806,] 
     182</stdout></cell><cell number="39"><input> 
     183print y[:36:6,1] 
     184</input><stdout> 
     185[ 0.355028, 0.339511, 0.324067, 0.308763, 0.293658, 0.278806,] 
     186</stdout></cell><cell number="40"><input> 
     187print y2[:36:6,1] 
     188</input><stdout> 
     189[ 0.355028, 0.339511, 0.324067, 0.308763, 0.293658, 0.278806,] 
     190</stdout></cell></ipython-log><sheet id="chapter04"> 
    2191<title>Integration (integrate)</title> 
    3 <para>The integrate sub-package provides several integration techniques 
    4 including an ordinary differential equation integrator. An overview of the 
    5 module is provided by the help command: 
    6 </para> 
    7 <programlisting>XXX: ipython 
    8 &gt;&gt;&gt; help(integrate) 
    9 Methods for Integrating Functions 
    10  
    11   odeint        -- Integrate ordinary differential equations. 
    12   quad          -- General purpose integration. 
    13   dblquad       -- General purpose double integration. 
    14   tplquad       -- General purpose triple integration. 
    15   gauss_quad    -- Integrate func(x) using Gaussian quadrature of order n. 
    16   gauss_quadtol -- Integrate with given tolerance using Gaussian quadrature. 
    17  
    18   See the orthogonal module (integrate.orthogonal) for Gaussian 
    19      quadrature roots and weights. 
    20 </programlisting> 
     192<para>The integrate sub-package provides several integration techniques including an ordinary differential equation integrator. An overview of the module is provided by the help command: </para><ipython-block logid="default-log"><ipython-cell type="input" number="1"/><ipython-cell type="input" number="2"/><ipython-cell type="stdout" number="2"/></ipython-block> 
    21193<section> 
    22194<title>General integration (integrate.quad)</title> 
    …… 
    30202This could be computed using quad: 
    31203</para> 
    32 <programlisting>XXX: ipython 
    33 &gt;&gt;&gt; result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5) 
    34 &gt;&gt;&gt; print result 
    35 (1.1178179380783249, 7.8663172481899801e-09) 
    36  
    37 &gt;&gt;&gt; I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5)-4.0/27*sqrt(2)*sin(4.5)+ 
    38     sqrt(2*pi)*special.fresnl(3/sqrt(pi))[0]) 
    39 &gt;&gt;&gt; print I 
    40 1.117817938088701 
    41  
    42 &gt;&gt;&gt; print abs(result[0]-I) 
    43 1.03761443881e-11  
    44 </programlisting> 
     204<ipython-block logid="default-log"><ipython-cell type="input" number="3"/><ipython-cell type="input" number="4"/><ipython-cell type="stdout" number="4"/><ipython-cell type="input" number="5"/><ipython-cell type="input" number="6"/><ipython-cell type="stdout" number="6"/><ipython-cell type="input" number="7"/><ipython-cell type="stdout" number="7"/></ipython-block> 
    45205<para>The first argument to quad is a "callable" Python object (i.e a function, 
    46206method, or class instance). Notice the use of a lambda-function in this case as 
    …… 
    72232defining a new function vec_expint based on the routine quad:  
    73233</para> 
    74 <programlisting> XXX: ipython 
    75 &gt;&gt;&gt; from integrate import quad, Inf 
    76 &gt;&gt;&gt; def integrand(t,n,x): 
    77         return exp(-x*t) / t**n 
    78  
    79 &gt;&gt;&gt; def expint(n,x):  
    80         return quad(integrand, 1, Inf, args=(n, x))[0] 
    81  
    82 &gt;&gt;&gt; vec_expint = vectorize(expint) 
    83  
    84 &gt;&gt;&gt; vec_expint(3,arange(1.0,4.0,0.5)) 
    85 array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049]) 
    86 &gt;&gt;&gt; special.expn(3,arange(1.0,4.0,0.5)) 
    87 array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049]) 
    88 </programlisting> 
     234 
     235<ipython-block logid="default-log"><ipython-cell type="input" number="8"/><ipython-cell type="input" number="9"/><ipython-cell type="input" number="10"/><ipython-cell type="input" number="11"/><ipython-cell type="input" number="12"/><ipython-cell type="output" number="12"/><ipython-cell type="input" number="13"/><ipython-cell type="output" number="13"/></ipython-block> 
    89236<para>The function which is integrated can even use the quad argument (though 
    90237the error bound may underestimate the error due to possible numerical error in 
    …… 
    94241dx=\frac{1}{n}</ipython-equation>. 
    95242</para> 
    96 <programlisting>XXX: ipython 
    97 &gt;&gt;&gt; result = quad(lambda x: expint(3, x), 0, Inf) 
    98 &gt;&gt;&gt; print result 
    99 (0.33333333324560266, 2.8548934485373678e-09)   
    100  
    101 &gt;&gt;&gt; I3 = 1.0/3.0 
    102 &gt;&gt;&gt; print I3 
    103 0.333333333333    
    104  
    105 &gt;&gt;&gt; print I3 - result[0] 
    106 8.77306560731e-11      
    107 </programlisting> 
     243<ipython-block logid="default-log"><ipython-cell type="input" number="14"/><ipython-cell type="input" number="15"/><ipython-cell type="stdout" number="15"/><ipython-cell type="input" number="16"/><ipython-cell type="input" number="17"/><ipython-cell type="stdout" number="17"/><ipython-cell type="input" number="18"/><ipython-cell type="stdout" number="18"/></ipython-block> 
     244 
    108245<para>This last example shows that multiple integration can be handled using 
    109246repeated calls to quad. The mechanics of this for double and triple integration 
    …… 
    115252<ipython-inlineequation>I_{n}</ipython-inlineequation> is shown below: 
    116253</para> 
    117 <programlisting>XXX: ipython 
    118 &gt;&gt;&gt; from __future__ import nested_scopes 
    119 &gt;&gt;&gt; from integrate import quad, dblquad, Inf 
    120 &gt;&gt;&gt; def I(n): 
    121     return dblquad(lambda t, x: exp(-x*t)/t**n, 0, Inf, lambda x: 1, lambda x: Inf)  
    122  
    123 &gt;&gt;&gt; print I(4) 
    124 (0.25000000000435768, 1.0518245707751597e-09) 
    125 &gt;&gt;&gt; print I(3) 
    126 (0.33333333325010883, 2.8604069919261191e-09)  
    127 &gt;&gt;&gt; print I(2) 
    128 (0.49999999999857514, 1.8855523253868967e-09) 
    129 </programlisting> 
     254 
     255<ipython-block logid="default-log"><ipython-cell type="input" number="19"/><ipython-cell type="input" number="20"/><ipython-cell type="input" number="21"/><ipython-cell type="stdout" number="21"/><ipython-cell type="input" number="22"/><ipython-cell type="stdout" number="22"/><ipython-cell type="input" number="23"/><ipython-cell type="stdout" number="23"/></ipython-block> 
    130256</section> 
    131257<section> 
    132258<title>Gaussian quadrature (integrate.gauss_quadtol)</title> 
    133 <para>A few functions are also provided in order to perform simple Gaussian 
    134 quadrature over a fixed interval. The first is fixed_quad which performs 
    135 fixed-order Gaussian quadrature. The second function is quadrature which 
    136 performs Gaussian quadrature of multiple orders until the difference in the 
    137 integral estimate is beneath some tolerance supplied by the user. These 
    138 functions both use the module special.orthogonal which can calculate the roots 
    139 and quadrature weights of a large variety of orthogonal polynomials (the 
    140 polynomials themselves are available as special functions returning instances of 
    141 the polynomial class --- e.g. special.legendre).  
    142 </para> 
    143 </section> 
     259<para>A few functions are also provided in order to perform simple Gaussian quadrature over a fixed interval. The first is fixed_quad which performs fixed-order Gaussian quadrature. The second function is quadrature which performs Gaussian quadrature of multiple orders until the difference in the integral estimate is beneath some tolerance supplied by the user. These functions both use the module special.orthogonal which can calculate the roots and quadrature weights of a large variety of orthogonal polynomials (the polynomials themselves are available as special functions returning instances of the polynomial class --- e.g. special.legendre). </para></section> 
    144260<section> 
    145261<title>Integrating using samples</title> 
    146 <para>There are three functions for computing integrals given only samples: 
    147 trapz, simps, and romb. The first two functions use Newton-Coates formulas of 
    148 order 1 and 2 respectively to perform integration. These two functions can 
    149 handle, non-equally-spaced samples. The trapezoidal rule approximates the 
    150 function as a straight line between adjacent points, while Simpson's rule 
    151 approximates the function between three adjacent points as a parabola.  
    152 </para> 
    153 <para>If the samples are equally-spaced and the number of samples available is 
     262<para>There are three functions for computing integrals given only samples: trapz, simps, and romb. The first two functions use Newton-Coates formulas of order 1 and 2 respectively to perform integration. These two functions can handle, non-equally-spaced samples. The trapezoidal rule approximates the function as a straight line between adjacent points, while Simpson's rule approximates the function between three adjacent points as a parabola. </para><para>If the samples are equally-spaced and the number of samples available is 
    154263<ipython-inlineequation>2^{k}+1</ipython-inlineequation> for some integer  
    155264<ipython-inlineequation>k</ipython-inlineequation>, then Romberg integration can be used to obtain 
    …… 
    223332<ipython-inlineequation>\mathbf{A}\left(t\right)</ipython-inlineequation> and its integral do not commute. 
    224333</para> 
    225 <para> 
    226 There are many optional inputs and outputs available when using odeint which can 
    227 help tune the solver. These additional inputs and outputs are not needed much of 
    228 the time, however, and the three required input arguments and the output 
    229 solution suffice. The required inputs are the function defining the derivative, 
    230 fprime, the initial conditions vector, y0, and the time points to obtain a 
    231 solution, t, (with the initial value point as the first element of this 
    232 sequence). The output to odeint is a matrix where each row contains the solution 
    233 vector at each requested time point (thus, the initial conditions are given in 
    234 the first output row).  
    235 </para> 
    236 <para> 
     334<para>There are many optional inputs and outputs available when using odeint which can help tune the solver. These additional inputs and outputs are not needed much of the time, however, and the three required input arguments and the output solution suffice. The required inputs are the function defining the derivative, fprime, the initial conditions vector, y0, and the time points to obtain a solution, t, (with the initial value point as the first element of this sequence). The output to odeint is a matrix where each row contains the solution vector at each requested time point (thus, the initial conditions are given in the first output row). </para><para> 
    237335The following example illustrates the use of odeint including the usage of the 
    238336Dfun option which allows the user to specify a gradient (with respect to 
    …… 
    240338<ipython-inlineequation>\mathbf{f}\left(\mathbf{y},t\right)</ipython-inlineequation>. 
    241339</para> 
    242 <programlisting>XXX: ipython 
    243 &gt;&gt;&gt; from integrate import odeint 
    244 &gt;&gt;&gt; from special import gamma, airy 
    245 &gt;&gt;&gt; y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0) 
    246 &gt;&gt;&gt; y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0) 
    247 &gt;&gt;&gt; y0 = [y0_0, y1_0] 
    248 &gt;&gt;&gt; def func(y, t): 
    249         return [t*y[1],y[0]] 
    250  
    251 &gt;&gt;&gt; def gradient(y,t): 
    252         return [[0,t],[1,0]] 
    253  
    254 &gt;&gt;&gt; x = arange(0,4.0, 0.01) 
    255 &gt;&gt;&gt; t = x 
    256 &gt;&gt;&gt; ychk = airy(x)[0] 
    257 &gt;&gt;&gt; y = odeint(func, y0, t) 
    258 &gt;&gt;&gt; y2 = odeint(func, y0, t, Dfun=gradient) 
    259  
    260 &gt;&gt;&gt; import sys 
    261 &gt;&gt;&gt; sys.float_output_precision = 6 
    262 &gt;&gt;&gt; print ychk[:36:6] 
    263 [ 0.355028  0.339511  0.324068  0.308763  0.293658  0.278806] 
    264  
    265 &gt;&gt;&gt; print y[:36:6,1] 
    266 [ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806] 
    267  
    268 &gt;&gt;&gt; print y2[:36:6,1] 
    269 [ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806] 
    270 </programlisting> 
     340<ipython-block logid="default-log"><ipython-cell type="input" number="24"/><ipython-cell type="input" number="25"/><ipython-cell type="input" number="26"/><ipython-cell type="input" number="27"/><ipython-cell type="input" number="28"/><ipython-cell type="input" number="29"/><ipython-cell type="input" number="30"/><ipython-cell type="input" number="31"/><ipython-cell type="input" number="32"/><ipython-cell type="input" number="33"/><ipython-cell type="input" number="34"/><ipython-cell type="input" number="35"/><ipython-cell type="input" number="36"/><ipython-cell type="input" number="37"/><ipython-cell type="input" number="38"/><ipython-cell type="stdout" number="38"/><ipython-cell type="input" number="39"/><ipython-cell type="stdout" number="39"/><ipython-cell type="input" number="40"/><ipython-cell type="stdout" number="40"/></ipython-block> 
     341 
    271342 
    272343</section> 

  • nbdoc/trunk/scipy_tutorial/chapter05.pybk

    r901 r902  
    1 <notebook version="1"><head/><ipython-log id="default-log"/><sheet id="chapter05"> 
     1<notebook version="1"><head/><ipython-log id="default-log"><cell number="0"><input>   
     2############DO NOT EDIT THIS CELL############   
     3from pylab import *   
     4switch_backend('WXAgg')   
     5ion()   
     6</input></cell><cell number="1"><input> 
     7from scipy.optimize import fmin 
     8</input></cell><cell number="2"><input> 
     9import scipy; from scipy import * 
     10 
     11</input></cell><cell number="3"><input> 
     12def rosen(x): 
     13  return sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) 
     14</input></cell><cell number="4"><input> 
     15x0 = [1.3, 0.7, 0.8, 1.9, 1.2] 
     16</input></cell><cell number="5"><input> 
     17xopt = fmin(rosen, x0) 
     18</input><stdout> 
     19Optimization terminated successfully. 
     20         Current function value: 0.000066 
     21         Iterations: 141 
     22         Function evaluations: 243 
     23</stdout></cell><cell number="6"><input> 
     24print xopt 
     25</input><stdout> 
     26[ 0.99910115, 0.99820923, 0.99646346, 0.99297555, 0.98600385,] 
     27</stdout></cell><cell number="7"><input> 
     28rosen(ones(5)) 
     29</input><output> 
     300.0 
     31</output></cell><cell number="8"><input> 
     32def rosen_der(x): 
     33  xm = x[1:-1] 
     34  xm_m1 = x[:-2] 
     35  xm_p1 = x[2:] 
     36  der = zeros(x.shape, x.typecode()) 
     37  der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm) 
     38  der[0] = -400*x[0]*(x[1] - x[0]**2) - 2*(1-x[0]) 
     39  der[-1] = 200*(x[-1]-x[-2]**2) 
     40  return der 
     41   </