Ticket #12922: 12922.patch

File 12922.patch, 14.6 KB (added by mmarco, 8 years ago)

Newer version, solves the pynac registry pollution without manipulating it

  • sage/calculus/functional.py

    # HG changeset patch
    # User Miguel Marco <mmarco@unizar.es>
    # Date 1370626733 -7200
    # Node ID 2d0de9bd9405fa847b097af5835eb11c5704635e
    # Parent  4536d5b2236a87d16e81e948c74f16ce8d8ef715
    Trac #12922 implicit derivative
    
    diff --git a/sage/calculus/functional.py b/sage/calculus/functional.py
    a b  
    457457        Traceback (most recent call last):
    458458        ...
    459459        ValueError: Expression F contains no y terms
     460       
     461   
     462    TESTS::
     463   
     464        sage: psr=copy(sage.symbolic.ring.pynac_symbol_registry)
     465        sage: var('x,y')
     466        (x, y)
     467        sage: psr=copy(sage.symbolic.ring.pynac_symbol_registry)
     468        sage: implicit_derivative(y,x,x^6*y^5,3)
     469        -792/125*y/x^3 + 12/25*(15*x^4*y^5 + 28*x^3*y^5)/(x^6*y^4) - 36/125*(20*x^5*y^4 + 43*x^4*y^4)/(x^7*y^3)
     470        sage: psr==sage.symbolic.ring.pynac_symbol_registry
     471        True
     472
     473       
    460474    """
    461475
    462476    if not F.has(Y):
    463477        raise ValueError, "Expression F contains no %s terms" % str(Y)
     478    import copy
     479    import sage
     480    psr=copy.copy(sage.symbolic.ring.pynac_symbol_registry)
    464481    x=SR.symbol('x')
    465482    yy=SR.symbol('yy')
    466483    y=SymbolicFunction('y',1)(x)
     
    479496        for j in range(n+1-i):
    480497            di[f.diff(x,i,yy,j)(x=x,yy=y)]=F.diff(X,i,Y,j)
    481498            S=S.subs_expr(di)
     499    sage.symbolic.ring.pynac_symbol_registry=psr
    482500    return S
  • sage/calculus/all.py

    # HG changeset patch
    # User Miguel Marco <mmarco@unizar.es>
    # Date 1370707293 -7200
    # Node ID 1b144e1544964b893bf7446620f7c838038e1ea9
    # Parent  01cb5f3a60abf5c1aa764987eb2fbd08ccd4e151
    Trac #12922 implicit_derivative
    
    diff --git a/sage/calculus/all.py b/sage/calculus/all.py
    a b  
    44
    55from functional import (diff, derivative,
    66                        expand,
    7                         taylor, simplify)
     7                        taylor, simplify, implicit_derivative)
    88
    99from functions import (wronskian,jacobian)
    1010
    1111from desolvers import (desolve, desolve_laplace, desolve_system,
    12                        eulers_method, eulers_method_2x2, 
    13                        eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4, 
     12                       eulers_method, eulers_method_2x2,
     13                       eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4,
    1414                       desolve_odeint)
    1515
    1616from var import (var, function, clear_vars)
     
    3535    - a symbolic expression.
    3636
    3737    EXAMPLES::
    38    
     38
    3939        sage: a = symbolic_expression(3/2); a
    4040        3/2
    4141        sage: type(a)
     
    6161        sage: symbolic_expression(E)
    6262        x*y + y^2 + y == x^3 + x^2 - 10*x - 10
    6363        sage: symbolic_expression(E) in SR
    64         True         
    65        
     64        True
     65
    6666    If x is a list or tuple, create a vector of symbolic expressions::
    67    
     67
    6868        sage: v=symbolic_expression([x,1]); v
    6969        (x, 1)
    7070        sage: v.base_ring()
    7171        Symbolic Ring
    7272        sage: v=symbolic_expression((x,1)); v
    7373        (x, 1)
    74         sage: v.base_ring()               
     74        sage: v.base_ring()
    7575        Symbolic Ring
    7676        sage: v=symbolic_expression((3,1)); v
    7777        (3, 1)
    78         sage: v.base_ring()               
     78        sage: v.base_ring()
    7979        Symbolic Ring
    8080        sage: E = EllipticCurve('15a'); E
    8181        Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 10*x - 10 over Rational Field
  • sage/calculus/functional.py

    diff --git a/sage/calculus/functional.py b/sage/calculus/functional.py
    a b  
    3232
    3333from calculus import SR
    3434from sage.symbolic.expression import Expression
     35from sage.symbolic.function_factory import SymbolicFunction
    3536
    3637def simplify(f):
    3738    r"""
    3839    Simplify the expression `f`.
    39    
     40
    4041    EXAMPLES: We simplify the expression `i + x - x`.
    41    
     42
    4243    ::
    43    
     44
    4445        sage: f = I + x - x; simplify(f)
    4546        I
    46    
     47
    4748    In fact, printing `f` yields the same thing - i.e., the
    4849    simplified form.
    4950    """
     
    5556def derivative(f, *args, **kwds):
    5657    """
    5758    The derivative of `f`.
    58    
     59
    5960    Repeated differentiation is supported by the syntax given in the
    6061    examples below.
    61    
     62
    6263    ALIAS: diff
    63    
     64
    6465    EXAMPLES: We differentiate a callable symbolic function::
    65    
     66
    6667        sage: f(x,y) = x*y + sin(x^2) + e^(-x)
    67         sage: f 
     68        sage: f
    6869        (x, y) |--> x*y + e^(-x) + sin(x^2)
    6970        sage: derivative(f, x)
    7071        (x, y) |--> 2*x*cos(x^2) + y - e^(-x)
    7172        sage: derivative(f, y)
    7273        (x, y) |--> x
    73    
     74
    7475    We differentiate a polynomial::
    75    
     76
    7677        sage: t = polygen(QQ, 't')
    7778        sage: f = (1-t)^5; f
    7879        -t^5 + 5*t^4 - 10*t^3 + 10*t^2 - 5*t + 1
     
    8687        -20*t^3 + 60*t^2 - 60*t + 20
    8788        sage: derivative(f, 2)
    8889        -20*t^3 + 60*t^2 - 60*t + 20
    89    
     90
    9091    We differentiate a symbolic expression::
    91    
     92
    9293        sage: var('a x')
    9394        (a, x)
    9495        sage: f = exp(sin(a - x^2))/x
     
    9697        -2*e^(sin(-x^2 + a))*cos(-x^2 + a) - e^(sin(-x^2 + a))/x^2
    9798        sage: derivative(f, a)
    9899        e^(sin(-x^2 + a))*cos(-x^2 + a)/x
    99    
     100
    100101    Syntax for repeated differentiation::
    101    
     102
    102103        sage: R.<u, v> = PolynomialRing(QQ)
    103104        sage: f = u^4*v^5
    104105        sage: derivative(f, u)
    105106        4*u^3*v^5
    106107        sage: f.derivative(u)   # can always use method notation too
    107108        4*u^3*v^5
    108    
     109
    109110    ::
    110    
     111
    111112        sage: derivative(f, u, u)
    112113        12*u^2*v^5
    113114        sage: derivative(f, u, u, u)
    114115        24*u*v^5
    115116        sage: derivative(f, u, 3)
    116117        24*u*v^5
    117    
     118
    118119    ::
    119    
     120
    120121        sage: derivative(f, u, v)
    121122        20*u^3*v^4
    122123        sage: derivative(f, u, 2, v)
     
    139140def integral(f, *args, **kwds):
    140141    r"""
    141142    The integral of `f`.
    142    
     143
    143144    EXAMPLES::
    144    
     145
    145146        sage: integral(sin(x), x)
    146147        -cos(x)
    147148        sage: integral(sin(x)^2, x, pi, 123*pi/2)
    148149        121/4*pi
    149150        sage: integral( sin(x), x, 0, pi)
    150151        2
    151    
     152
    152153    We integrate a symbolic function::
    153    
     154
    154155        sage: f(x,y,z) = x*y/z + sin(z)
    155156        sage: integral(f, z)
    156157        (x, y, z) |--> x*y*log(z) - cos(z)
    157    
     158
    158159    ::
    159    
     160
    160161        sage: var('a,b')
    161162        (a, b)
    162163        sage: assume(b-a>0)
    163164        sage: integral( sin(x), x, a, b)
    164165        cos(a) - cos(b)
    165166        sage: forget()
    166    
     167
    167168    ::
    168    
     169
    169170        sage: integral(x/(x^3-1), x)
    170171        1/3*sqrt(3)*arctan(1/3*(2*x + 1)*sqrt(3)) + 1/3*log(x - 1) - 1/6*log(x^2 + x + 1)
    171    
     172
    172173    ::
    173    
     174
    174175        sage: integral( exp(-x^2), x )
    175176        1/2*sqrt(pi)*erf(x)
    176    
     177
    177178    We define the Gaussian, plot and integrate it numerically and
    178179    symbolically::
    179    
     180
    180181        sage: f(x) = 1/(sqrt(2*pi)) * e^(-x^2/2)
    181182        sage: P = plot(f, -4, 4, hue=0.8, thickness=2)
    182183        sage: P.show(ymin=0, ymax=0.4)
     
    184185        (0.99993665751633376, 1.1101527003413533e-14)
    185186        sage: integrate(f, x)
    186187        x |--> 1/2*erf(1/2*sqrt(2)*x)
    187    
     188
    188189    You can have Sage calculate multiple integrals. For example,
    189190    consider the function `exp(y^2)` on the region between the
    190191    lines `x=y`, `x=1`, and `y=0`. We find the
    191192    value of the integral on this region using the command::
    192    
     193
    193194        sage: area = integral(integral(exp(y^2),x,0,y),y,0,1); area
    194195        1/2*e - 1/2
    195196        sage: float(area)
    196197        0.859140914229522...
    197    
     198
    198199    We compute the line integral of `\sin(x)` along the arc of
    199200    the curve `x=y^4` from `(1,-1)` to
    200201    `(1,1)`::
    201    
     202
    202203        sage: t = var('t')
    203204        sage: (x,y) = (t^4,t)
    204205        sage: (dx,dy) = (diff(x,t), diff(y,t))
    205206        sage: integral(sin(x)*dx, t,-1, 1)
    206207        0
    207208        sage: restore('x,y')   # restore the symbolic variables x and y
    208    
     209
    209210    Sage is unable to do anything with the following integral::
    210    
     211
    211212        sage: integral( exp(-x^2)*log(x), x )
    212213        integrate(e^(-x^2)*log(x), x)
    213    
     214
    214215    Note, however, that::
    215    
     216
    216217        sage: integral( exp(-x^2)*ln(x), x, 0, oo)
    217218        -1/4*(euler_gamma + 2*log(2))*sqrt(pi)
    218    
     219
    219220    This definite integral is easy::
    220    
     221
    221222        sage: integral( ln(x)/x, x, 1, 2)
    222223        1/2*log(2)^2
    223    
     224
    224225    Sage can't do this elliptic integral (yet)::
    225    
     226
    226227        sage: integral(1/sqrt(2*t^4 - 3*t^2 - 2), t, 2, 3)
    227228        integrate(1/sqrt(2*t^4 - 3*t^2 - 2), t, 2, 3)
    228    
     229
    229230    A double integral::
    230    
     231
    231232        sage: y = var('y')
    232233        sage: integral(integral(x*y^2, x, 0, y), y, -2, 2)
    233234        32/5
    234    
     235
    235236    This illustrates using assumptions::
    236    
     237
    237238        sage: integral(abs(x), x, 0, 5)
    238239        25/2
    239240        sage: a = var("a")
     
    251252        sage: integral(abs(x)*x, x, 0, a)
    252253        1/3*a^3
    253254        sage: forget()      # forget the assumptions.
    254    
     255
    255256    We integrate and differentiate a huge mess::
    256    
     257
    257258        sage: f = (x^2-1+3*(1+x^2)^(1/3))/(1+x^2)^(2/3)*x/(x^2+2)^2
    258259        sage: g = integral(f, x)
    259260        sage: h = f - diff(g, x)
    260    
     261
    261262    ::
    262    
     263
    263264        sage: [float(h(i)) for i in range(5)] #random
    264        
     265
    265266        [0.0,
    266267         -1.1102230246251565e-16,
    267268         -5.5511151231257827e-17,
     
    287288    r"""
    288289    Return the limit as the variable `v` approaches `a`
    289290    from the given direction.
    290    
     291
    291292    ::
    292    
     293
    293294                limit(expr, x = a)
    294295                limit(expr, x = a, dir='above')
    295                
    296    
     296
     297
    297298    INPUT:
    298    
     299
    299300    - ``dir`` - (default: None); dir may have the value
    300301       'plus' (or 'above') for a limit from above, 'minus' (or 'below')
    301302       for a limit from below, or may be omitted (implying a two-sided
    302303       limit is to be computed).
    303    
     304
    304305    - ``taylor`` - (default: False); if True, use Taylor
    305306       series, which allows more limits to be computed (but may also
    306307       crash in some obscure cases due to bugs in Maxima).
    307    
     308
    308309    - ``\*\*argv`` - 1 named parameter
    309    
     310
    310311    ALIAS: You can also use lim instead of limit.
    311    
     312
    312313    EXAMPLES::
    313    
     314
    314315        sage: limit(sin(x)/x, x=0)
    315316        1
    316317        sage: limit(exp(x), x=oo)
     
    323324        -1/2
    324325        sage: limit((tan(sin(x)) - sin(tan(x)))/x^7, taylor=True, x=0)
    325326        1/30
    326    
     327
    327328    Sage does not know how to do this limit (which is 0), so it returns
    328329    it unevaluated::
    329    
     330
    330331        sage: lim(exp(x^2)*(1-erf(x)), x=infinity)
    331332        -limit((erf(x) - 1)*e^(x^2), x, +Infinity)
    332333    """
     
    342343    variable `v` around the point `a`, containing terms
    343344    through `(x - a)^n`. Functions in more variables are also
    344345    supported.
    345    
     346
    346347    INPUT:
    347    
     348
    348349    - ``*args`` - the following notation is supported
    349    
     350
    350351    - ``x, a, n`` - variable, point, degree
    351    
     352
    352353    - ``(x, a), (y, b), ..., n`` - variables with points, degree of polynomial
    353    
     354
    354355    EXAMPLES::
    355    
     356
    356357        sage: var('x,k,n')
    357358        (x, k, n)
    358359        sage: taylor (sqrt (1 - k^2*sin(x)^2), x, 0, 6)
    359360        -1/720*(45*k^6 - 60*k^4 + 16*k^2)*x^6 - 1/24*(3*k^4 - 4*k^2)*x^4 - 1/2*k^2*x^2 + 1
    360        
     361
    361362    ::
    362363
    363364        sage: taylor ((x + 1)^n, x, 0, 4)
    364365        1/24*(n^4 - 6*n^3 + 11*n^2 - 6*n)*x^4 + 1/6*(n^3 - 3*n^2 + 2*n)*x^3 + 1/2*(n^2 - n)*x^2 + n*x + 1
    365        
     366
    366367    ::
    367368
    368369        sage: taylor ((x + 1)^n, x, 0, 4)
    369370        1/24*(n^4 - 6*n^3 + 11*n^2 - 6*n)*x^4 + 1/6*(n^3 - 3*n^2 + 2*n)*x^3 + 1/2*(n^2 - n)*x^2 + n*x + 1
    370        
     371
    371372    Taylor polynomial in two variables::
    372373
    373374        sage: x,y=var('x y'); taylor(x*y^3,(x,1),(y,-1),4)
    374         (y + 1)^3*(x - 1) + (y + 1)^3 - 3*(y + 1)^2*(x - 1) - 3*(y + 1)^2 + 3*(y + 1)*(x - 1) - x + 3*y + 3 
     375        (y + 1)^3*(x - 1) + (y + 1)^3 - 3*(y + 1)^2*(x - 1) - 3*(y + 1)^2 + 3*(y + 1)*(x - 1) - x + 3*y + 3
    375376    """
    376377    if not isinstance(f, Expression):
    377378        f = SR(f)
     
    380381def expand(x, *args, **kwds):
    381382    """
    382383    EXAMPLES::
    383    
     384
    384385        sage: a = (x-1)*(x^2 - 1); a
    385386        (x - 1)*(x^2 - 1)
    386387        sage: expand(a)
    387388        x^3 - x^2 - x + 1
    388    
     389
    389390    You can also use expand on polynomial, integer, and other
    390391    factorizations::
    391    
     392
    392393        sage: x = polygen(ZZ)
    393394        sage: F = factor(x^12 - 1); F
    394395        (x - 1) * (x + 1) * (x^2 - x + 1) * (x^2 + 1) * (x^2 + x + 1) * (x^4 - x^2 + 1)
     
    400401        3^2 * 223
    401402        sage: expand(F)
    402403        2007
    403    
     404
    404405    Note: If you want to compute the expanded form of a polynomial
    405406    arithmetic operation quickly and the coefficients of the polynomial
    406407    all lie in some ring, e.g., the integers, it is vastly faster to
    407408    create a polynomial ring and do the arithmetic there.
    408    
     409
    409410    ::
    410    
     411
    411412        sage: x = polygen(ZZ)      # polynomial over a given base ring.
    412413        sage: f = sum(x^n for n in range(5))
    413414        sage: f*f                  # much faster, even if the degree is huge
    414415        x^8 + 2*x^7 + 3*x^6 + 4*x^5 + 5*x^4 + 4*x^3 + 3*x^2 + 2*x + 1
    415    
     416
    416417    TESTS::
    417    
     418
    418419        sage: t1 = (sqrt(3)-3)*(sqrt(3)+1)/6;
    419420        sage: tt1 = -1/sqrt(3);
    420421        sage: t2 = sqrt(3)/6;
     
    433434        return x.expand(*args, **kwds)
    434435    except AttributeError:
    435436        return x
     437
     438
     439def implicit_derivative(Y,X,F,n=1):
     440    """
     441    Computes the n'th derivative of Y with respect to X given implicitly by F.
     442
     443    EXAMPLES:
     444
     445    ::
     446
     447        sage: var('x,y')
     448        (x, y)
     449        sage: implicit_derivative(y,x,cos(x)*sin(y))
     450        sin(x)*sin(y)/(cos(x)*cos(y))
     451        sage: implicit_derivative(y,x,x*y^2,3)
     452        -1/4*(y + 2*y/x)/x^2 + 1/4*(2*y^2/x - y^2/x^2)/(x*y) - 3/4*y/x^3
     453
     454    It is an error to not include a Y term in the expression F::
     455
     456        sage: implicit_derivative(y,x,cos(x)*sin(x))
     457        Traceback (most recent call last):
     458        ...
     459        ValueError: Expression F contains no y terms
     460
     461
     462    TESTS::
     463
     464        sage: psr=copy(sage.symbolic.ring.pynac_symbol_registry)
     465        sage: var('x,y')
     466        (x, y)
     467        sage: psr=copy(sage.symbolic.ring.pynac_symbol_registry)
     468        sage: implicit_derivative(y,x,x^6*y^5,3)
     469        -792/125*y/x^3 + 12/25*(15*x^4*y^5 + 28*x^3*y^5)/(x^6*y^4) - 36/125*(20*x^5*y^4 + 43*x^4*y^4)/(x^7*y^3)
     470        sage: psr==sage.symbolic.ring.pynac_symbol_registry
     471        True
     472
     473
     474    """
     475
     476    if not F.has(Y):
     477        raise ValueError, "Expression F contains no %s terms" % str(Y)
     478    x=SR.symbol()
     479    yy=SR.symbol()
     480    y=SymbolicFunction('y',1)(x)
     481    f=SymbolicFunction('f',2)(x,yy)
     482    Fx=f.diff(x)
     483    Fy=f.diff(yy)
     484    G=-(Fx/Fy)
     485    G=G.subs({yy:y})
     486    di={y.diff(x):-F.diff(X)/F.diff(Y)}
     487    R=G
     488    S=G.diff(x,n-1)
     489    for i in range(n+1):
     490        di[y.diff(x,i+1).subs({x:x})]=R
     491        S=S.subs(di)
     492        R=G.diff(x,i)
     493        for j in range(n+1-i):
     494            di[f.diff(x,i,yy,j).subs({x:x,yy:y})]=F.diff(X,i,Y,j)
     495            S=S.subs(di)
     496    return S