# Ticket #8132: trac_8132.patch

File trac_8132.patch, 82.0 KB (added by Robert Marik, 13 years ago)
• ## doc/en/constructions/calculus.rst

# HG changeset patch
# User Robert Marik <marik@mendelu.cz>
# Date 1264898007 -3600
# Node ID ea16778cc7ae89ebf12c73c371e2fdc0cbab046e
# Parent  b89cf68ccc1ca1cfc9f29e86e9b61770b41d7df7
Ticket #8132 - fix documentation related to ODE solvers

diff -r b89cf68ccc1c -r ea16778cc7ae doc/en/constructions/calculus.rst
 a sage: latex(f.diff(x)) k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + w x^{3} e^{\left(k x\right)} \cos\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right) If you type view(f.diff('x')) another window will open up If you type view(f.diff(x)) another window will open up displaying the compiled output. In the notebook, you can enter :: f = maxima('x^3 * %e^(k*x) * sin(w*x)') show(f) show(f.diff('x')) var('x k w') f = x^3 * e^(k*x) * sin(w*x) show(f) show(f.diff(x)) into a cell and press shift-enter for a similar result. You can also call Maxima indirectly using the commands also differentiate and integrate using the commands :: Ordinary differential equations =============================== Symbolically solving ODEs can be done using 's interface with Maxima. Numerical solution of ODEs can be done using 's interface Symbolically solving ODEs can be done using Sage interface with Maxima. See :: sage:desolvers? for available commands. Numerical solution of ODEs can be done using Sage interface with Octave (an experimental package), or routines in the GSL (Gnu Scientific Library). You can solve ODE's symbolically in Sage using the Maxima interface An example, how to solve ODE's symbolically in Sage using the Maxima interface (do not type the ...): :: sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1]) y=3*x-2*%e^(x-1) sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y']) y=%k1*%e^x+%k2*%e^-x+3*x sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y']) y=(%c-3*(-x-1)*%e^-x)*%e^x sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1]) y=-%e^-1*(5*%e^x-3*%e*x-3*%e) sage: y=function('y',x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1]) 3*x - 2*e^(x - 1) sage: desolve(diff(y,x,2) + 3*x == y, dvar = y) k1*e^x + k2*e^(-x) + 3*x sage: desolve(diff(y,x) + 3*x == y, dvar = y) (3*(x + 1)*e^(-x) + c)*e^x sage: desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand() 3*x - 5*e^(x - 1) + 3 sage: maxima.clear('x'); maxima.clear('f') sage: maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\ ...   ["x","f"], [0,1,2]) f(x)=x*%e^x+%e^x sage: f=function('f',x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2]) x*e^x + e^x sage: maxima.clear('x'); maxima.clear('f') sage: f = maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\ ...   ["x","f"]) sage: f f(x)=x*%e^x*('at('diff(f(x),x,1),x=0))-f(0)*x*%e^x+f(0)*%e^x sage: f.display2d() ! x  d        !                  x          x f(x) = x %e  (-- (f(x))!     ) - f(0) x %e  + f(0) %e dx       ! !x = 0 sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f) -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0) .. index: pair: differential equations; plot .. math:: x' = x+y, x(0) = 1; y' = x-y, y(0) = -1, for :math:0 <= t <= 2. Another way this system can be solved is to use the command desolve_system in calculus/examples. for :math:0 <= t <= 2. The same result can be obtained by using desolve_system_rk4:: sage: var('x y t') sage: P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2) sage: p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True) sage: p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red') sage: p1+p2 Another way this system can be solved is to use the command desolve_system. .. skip :: sage: attach os.environ['SAGE_ROOT'] + '/examples/calculus/desolvers.sage' sage: des = ["'diff(x(t),t)=x(t)+y(t)","'diff(y(t),t)=x(t)-y(t)"] sage: vars = ["t","x","y"] sage: ics = [0,1,-1] sage: desolve_system(des,vars,ics) [x(t)=cosh(2^(1/2)*t),y(t)=2*sinh(2^(1/2)*t)/sqrt(2)-cosh(2^(1/2)*t)] sage: t=var('t'); x=function('x',t); y=function('y',t) sage: des = [diff(x,t) == x+y, diff(y,t) == x-y] sage: desolve_system(des, [x,y], ics = [0, 1, -1]) [x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)] The output of this command is *not* a pair of functions.
• ## sage/calculus/calculus.py

diff -r b89cf68ccc1c -r ea16778cc7ae sage/calculus/calculus.py
 a If self has only one variable, then it returns the integral with respect to that variable. INPUT: If definite integration fails, it could be still possbile to evaluate definite integral using indefinite integral and Newton - Leibniz theorem (however, the user has to ensure that the indefinite integral is continuous on compact interval [a,b] and this theorem can be applied). INPUT: - v - (optional) a variable or variable name.  This can also be a tuple of the variable (optional) and endpoints (i.e., (x,0,1) or (0,1)). be a tuple of the variable (optional) and endpoints (i.e., (x,0,1) or (0,1)). - a - (optional) lower endpoint of definite integral - 'mathematica_free' - use http://integrals.wolfram.com/ EXAMPLES:: EXAMPLES:: sage: x = var('x') sage: h = sin(x)/(cos(x))^2
• ## sage/calculus/desolvers.py

diff -r b89cf68ccc1c -r ea16778cc7ae sage/calculus/desolvers.py
 a """ This file contains functions useful for solving differential equations which occur commonly in a 1st semester differential equations course. which occur commonly in a 1st semester differential equations course. For another numerical solver see :meth:ode_solver function and optional package Octave. * desolve -- Computes the "general solution" to a 1st or 2nd order ODE via Maxima. Commands: * desolve_laplace -- Solves an ODE using laplace transforms via Maxima. Initials conditions are optional. - desolve - Computes the "general solution" to a 1st or 2nd order ODE via Maxima. * desolve_system -- Solves any size system of 1st order odes using Maxima. Initials conditions are optional. - desolve_laplace - Solves an ODE using laplace transforms via Maxima. Initials conditions are optional. * desolve_rk4 -- Solves numerically IVP for one first order equation, returns list of points or plot - desolve_system - Solves any size system of 1st order odes using Maxima. Initials conditions are optional. * desolve_system_rk4 -- Solves numerically IVP for system of first order equations, returns list of points - desolve_rk4 - Solves numerically IVP for one first order equation, returns list of points or plot * eulers_method -- Approximate solution to a 1st order DE, presented as a table. - desolve_system_rk4 - Solves numerically IVP for system of first order equations, returns list of points * eulers_method_2x2 -- Approximate solution to a 1st order system of DEs, presented as a table. - eulers_method - Approximate solution to a 1st order DE, presented as a table. * eulers_method_2x2_plot -- Plots the sequence of points obtained from Euler's method. - eulers_method_2x2 - Approximate solution to a 1st order system of DEs, presented as a table. AUTHORS: David Joyner (3-2006)     -- Initial version of functions Marshall Hampton (7-2007) -- Creation of Python module and testing Robert Bradshaw (10-2008) -- Some interface cleanup. Robert Marik (10-2009) -- Some bugfixes and enhancements - eulers_method_2x2_plot - Plots the sequence of points obtained from Euler's method. AUTHORS: - David Joyner (3-2006) - Initial version of functions - Marshall Hampton (7-2007) - Creation of Python module and testing - Robert Bradshaw (10-2008) - Some interface cleanup. - Robert Marik (10-2009) - Some bugfixes and enhancements """ def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False): """ Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP. *Use* desolve?  *if the output in truncated in notebook.* INPUT: de    -- an expression or equation representing the ODE dvar  -- the dependent variable (hereafter called y) ics   -- (optional) the initial or boundary conditions for a first-order equation, specify the initial x and y for a second-order equation, specify the initial x, y, and dy/dx for a second-order boundary solution, specify initial and final x and y initial conditions gives error if the solution is not SymbolisEquation (as happens for example for clairot equation) ivar  -- (optional) the independent variable  (hereafter called x), - de - an expression or equation representing the ODE - dvar - the dependent variable (hereafter called y) - ics - (optional) the initial or boundary conditions - for a first-order equation, specify the initial x and y - for a second-order equation, specify the initial x, y, and dy/dx - for a second-order boundary solution, specify initial and final x and y initial conditions gives error if the solution is not SymbolisEquation (as happens for example for clairot equation) - ivar - (optional) the independent variable  (hereafter called x), which must be specified if there is more than one independent variable in the equation. show_method  -- (optional) if true, then Sage returns pair [solution, method], where method is the string describing method which has been used to get solution (Maxima uses the following order for first order equations: linear, separable, exact (including exact with integrating factor), homogeneous, bernoulli, generalized homogeneous) - use carefully in class, see below for the example of the equation which is separable but this property is not recognized by Maxima and equation is solved as exact. contrib_ode  -- (optional) if true, desolve allows to solve clairot, lagrange, riccati and some other equations. May take a long time and thus turned off by default. Initial conditions can be used only if the result is one SymbolicEquation (does not contain singular solution, for example) - show_method - (optional) if true, then Sage returns pair [solution, method], where method is the string describing method which has been used to get solution (Maxima uses the following order for first order equations: linear, separable, exact (including exact with integrating factor), homogeneous, bernoulli, generalized homogeneous) - use carefully in class, see below for the example of the equation which is separable but this property is not recognized by Maxima and equation is solved as exact. - contrib_ode - (optional) if true, desolve allows to solve clairot, lagrange, riccati and some other equations. May take a long time and thus turned off by default.  Initial conditions can be used only if the result is one SymbolicEquation (does not contain singular solution, for example) OUTPUT: In most cases returns SymbolicEquation which defines the solution implicitly. If the result is in the form y(x)=... (happens for linear eqs.), returns the right-hand side only. The possible constant solutions of separable ODE's are omitted. In most cases returns SymbolicEquation which defines the solution implicitly.  If the result is in the form y(x)=... (happens for linear eqs.), returns the right-hand side only.  The possible constant solutions of separable ODE's are omitted. EXAMPLES: EXAMPLES:: sage: x = var('x') sage: y = function('y', x) sage: desolve(diff(y,x) + y - 1, y) (c + e^x)*e^(-x) :: sage: f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f (e^10 + e^x)*e^(-x) :: sage: plot(f) We can also solve second-order differential equations. We can also solve second-order differential equations.:: sage: x = var('x') sage: y = function('y', x) sage: de = diff(y,x,2) - y == x sage: desolve(de, y) k1*e^x + k2*e^(-x) - x :: sage: f = desolve(de, y, [10,2,1]); f -x + 5*e^(-x + 10) + 7*e^(x - 10) sage: f(x=10) sage: diff(f,x)(x=10) 1 :: sage: de = diff(y,x,2) + y == 0 sage: desolve(de, y) k1*sin(x) + k2*cos(x) sage: desolve(de, y, [0,1,pi/2,4]) 4*sin(x) + cos(x) sage: desolve(y*diff(y,x)+sin(x)==0,y) -1/2*y(x)^2 == c - cos(x) Clairot equation: general and sungular solutions Clairot equation: general and singular solutions:: sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True) [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault'] For equations involving more variables we specify independent variable For equations involving more variables we specify independent variable:: sage: a,b,c,n=var('a b c n') sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True) [[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati'] Higher orded, not involving independent variable Higher orded, not involving independent variable:: sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand() 1/6*y(x)^3 + k1*y(x) == k2 + x sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True) [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx'] Separable equations - Sage returns solution in implicit form Separable equations - Sage returns solution in implicit form:: sage: desolve(diff(y,x)*sin(y) == cos(x),y) -cos(y(x)) == c + sin(x) sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1]) -cos(y(x)) == sin(x) - cos(1) - 1 Linear equation - Sage returns the expression on the right hand side only Linear equation - Sage returns the expression on the right hand side only:: sage: desolve(diff(y,x)+(y) == cos(x),y) 1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x) sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1]) 1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x) This ODE with separated variables is solved as exact Explanation - factor does not split exp(x-y) in maxima into exp(x)*exp(-y) This ODE with separated variables is solved as exact. Explanation - factor does not split e^{x-y} in Maxima into e^{x}e^{y}:: sage: desolve(diff(y,x)==exp(x-y),y,show_method=True) [-e^x + e^y(x) == c, 'exact'] You can solve Bessel equations. You can also use initial conditions, but you cannot put (sometimes desired) initial condition at x=0, since this point is singlar point of the eruation. Anyway, if the solution should be bounded at x=0, then k2=0. equation. Anyway, if the solution should be bounded at x=0, then k2=0.:: sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y) k1*bessel_j(2, x) + k2*bessel_y(2, x) This example asks (silly) question, whether sqrt(3) is an integer. has been reported to Maxima developers. # sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-3)*y==0,y) Difficult ODE produces error:: sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) # not tested Traceback (click to the left for traceback) ... NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." Difficult ODE produces error - moreover, takes a long time :: Produces error (difficult ODE) #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) #Traceback (click to the left for traceback) #... #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested #Produces error (difficult ODE) - moreover, takes a long time #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) Some more types od ODE's Some more types od ODE's:: sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True) [[y(x) == c + log(x), y(x) == c*e^x], 'factor'] :: sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True) [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange'] These two examples produce error (as expected, Maxima 5.18 cannot solve equations from initial conditions) Current Maxima 5.18 returns false answer in this case! These two examples produce error (as expected, Maxima 5.18 cannot solve equations from initial conditions). Current Maxima 5.18 returns false answer in this case!:: #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand() #Traceback (click to the left for traceback) #... #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) #Traceback (click to the left for traceback) #... #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand() # not tested Traceback (click to the left for traceback) ... NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." :: sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested Traceback (click to the left for traceback) ... NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." Second order linear ODE Second order linear ODE:: sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y) (k2*x + k1)*e^(-x) + 1/2*sin(x) 3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True) [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y) (k2*x + k1)*e^(-x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True) [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff'] This equation can be solved within Maxima but not within Sage needs assumptions assume(x>0,y>0) and works in maxima, but not in Sage # sage: assume(x>0) # sage: assume(y>0) # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True) This equation can be solved within Maxima but not within Sage. It needs assumptions assume(x>0,y>0) and works in Maxima, but not in Sage:: This equation can be solved within Maxima but not within Sage Perhaps upgrading Maxima help # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True) sage: assume(x>0) # not tested sage: assume(y>0) # not tested sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True) # not tested TESTS: Trac #6479 fixed Trac #6479 fixed:: sage: x = var('x') sage: y = function('y', x) sage: desolve( diff(y,x,x) == 0, y, [0,0,1]) x :: sage: desolve( diff(y,x,x) == 0, y, [0,1,1]) x + 1 AUTHOR: David Joyner (1-2006) Robert Bradshaw (10-2008) Robert Marik (10-2009) AUTHORS: - David Joyner (1-2006) - Robert Bradshaw (10-2008) - Robert Marik (10-2009) """ if is_SymbolicEquation(de): Solves an ODE using laplace transforms. Initials conditions are optional. INPUT: de    -- a lambda expression representing the ODE (eg, de = diff(y,x,2) == diff(y,x)+sin(x)) dvar  -- the dependent variable (eg y) ivar  -- (optional) the independent variable  (hereafter called x), which must be specified if there is more than one independent variable in the equation. ics   -- a list of numbers representing initial conditions, (eg, f(0)=1, f'(0)=2 is ics = [0,1,2]) - de - a lambda expression representing the ODE (eg, de = diff(y,x,2) == diff(y,x)+sin(x)) - dvar - the dependent variable (eg y) - ivar - (optional) the independent variable (hereafter called x), which must be specified if there is more than one independent variable in the equation. - ics - a list of numbers representing initial conditions, (eg, f(0)=1, f'(0)=2 is ics = [0,1,2]) OUTPUT: Solution of the ODE as symbolic expression EXAMPLES: Solution of the ODE as symbolic expression EXAMPLES:: sage: u=function('u',x) sage: eq = diff(u,x) - exp(-x) - u == 0 sage: desolve_laplace(eq,u) 1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x) We can use initial conditions We can use initial conditions:: sage: desolve_laplace(eq,u,ics=[0,3]) -1/2*e^(-x) + 7/2*e^x The initial conditions do not persist in the system (as they persisted in previous versions) in previous versions):: sage: desolve_laplace(eq,u) 1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x) :: sage: f=function('f', x) sage: eq = diff(f,x) + f == 0 sage: desolve_laplace(eq,f,[0,1]) e^(-x) :: sage: x = var('x') sage: f = function('f', x) sage: de = diff(f,x,x) - 2*diff(f,x) + f x*e^x + e^x TESTS: #4839 fixed Trac #4839 fixed:: sage: t=var('t') sage: x=function('x', t) e^(-t) + 1 sage: soln(t=3) e^(-3) + 1 AUTHORS: AUTHOR: David Joyner (1-2006,8-2007) Robert Marik (10-2009) - David Joyner (1-2006,8-2007) - Robert Marik (10-2009) """ #This is the original code from David Joyner (inputs and outputs strings) #maxima("de:"+de._repr_()+"=0;") Solves any size system of 1st order ODE's. Initials conditions are optional. INPUT: des   -- list of ODEs vars  -- list of dependent variables ics   -- (optional) list of initial values for ivar and vars ivar  -- (optional) the independent variable, which must be specified if there is more than one independent variable in the equation. - des - list of ODEs - vars - list of dependent variables - ics - (optional) list of initial values for ivar and vars - ivar - (optional) the independent variable, which must be specified if there is more than one independent variable in the equation. EXAMPLES: EXAMPLES:: sage: t = var('t') sage: x = function('x', t) sage: y = function('y', t) [x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1, y(t) == (x(0) - 1)*sin(t) + (y(0) - 1)*cos(t) + 1] Now we give some initial conditions: Now we give some initial conditions:: sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol [x(t) == -sin(t) + 1, y(t) == cos(t) + 1] sage: solnx, solny = sol[0].rhs(), sol[1].rhs() sage: plot([solnx,solny],(0,1)) sage: parametric_plot((solnx,solny),(0,1)) AUTHOR: Robert Bradshaw (10-2008) AUTHORS: - Robert Bradshaw (10-2008) """ ivars = set([]) for i, de in enumerate(des): Now type show(P1), show(P2) to view these. AUTHOR: David Joyner (3-2006, 8-2007) AUTHORS: David Joyner (3-2006, 8-2007) """ d = len(des) dess = [de._maxima_init_() + "=0" for de in des] def eulers_method(f,x0,y0,h,x1,method="table"): """ This implements Euler's method for finding numerically the solution of the 1st order ODE y' = f(x,y), y(a)=c. The "x" column of the table increments from x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the This implements Euler's method for finding numerically the solution of the 1st order ODE y' = f(x,y), y(a)=c. The "x" column of the table increments from x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the last column. *For pedagogical purposes only.* EXAMPLES: EXAMPLES:: sage: from sage.calculus.desolvers import eulers_method sage: x,y = PolynomialRing(QQ,2,"xy").gens() sage: eulers_method(5*x+y-5,0,1,1/2,1) 0                    1                   -2 1/2                   -1                 -7/4 1                -11/4                -11/8 :: sage: x,y = PolynomialRing(QQ,2,"xy").gens() sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none") [[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]] :: sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: x,y = PolynomialRing(RR,2,"xy").gens() sage: eulers_method(5*x+y-5,0,1,1/2,1,method="None") [[0, 1], [1/2, -1.0], [1, -2.7], [3/2, -4.0]] :: sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: x,y=PolynomialRing(RR,2,"xy").gens() sage: eulers_method(5*x+y-5,0,1,1/2,1) 0                    1                 -2.0 1/2                 -1.0                 -1.7 1                 -2.7                 -1.3 :: sage: x,y=PolynomialRing(QQ,2,"xy").gens() sage: eulers_method(5*x+y-5,1,1,1/3,2) x                    y                  h*f(x,y) 4/3                  4/3                    1 5/3                  7/3                 17/9 2                 38/9                83/27 :: sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none") [[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]] :: sage: pts = eulers_method(5*x+y-5,0,1,1/2,1,method="none") sage: P1 = list_plot(pts) sage: P2 = line(pts) sage: (P1+P2).show() AUTHOR: David Joyner AUTHORS: - David Joyner """ if method=="table": print "%10s %20s %25s"%("x","y","h*f(x,y)") def eulers_method_2x2(f,g, t0, x0, y0, h, t1,method="table"): """ This implements Euler's method for finding numerically the solution of the 1st order system of two ODEs x' = f(t, x, y), x(t0)=x0. y' = g(t, x, y), y(t0)=y0. The "t" column of the table increments from t0 to t1 by h (so (t1-t0)/h must be an integer). In the "x" column, the new x-value equals the old x-value plus the corresponding entry in the next (third) column. In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the next (last) column. This implements Euler's method for finding numerically the solution of the 1st order system of two ODEs x' = f(t, x, y), x(t0)=x0. y' = g(t, x, y), y(t0)=y0. The "t" column of the table increments from t_0 to t_1 by h (so \\frac{t_1-t_0}{h} must be an integer). In the "x" column, the new x-value equals the old x-value plus the corresponding entry in the next (third) column.  In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the next (last) column. *For pedagogical purposes only.* EXAMPLES: EXAMPLES:: sage: from sage.calculus.desolvers import eulers_method_2x2 sage: t, x, y = PolynomialRing(QQ,3,"txy").gens() sage: f = x+y+t; g = x-y sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,method="none") [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]] :: sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1) t                    x                h*f(t,x,y)                    y           h*g(t,x,y) 0                    0                         0                    0                    0 1/3                    0                       1/9                    0                    0 2/3                  1/9                      7/27                    0                 1/27 1                10/27                     38/81                 1/27                  1/9 :: sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: t,x,y=PolynomialRing(RR,3,"txy").gens() sage: f = x+y+t; g = x-y 2/3                 0.13                      0.29                 0.00                0.043 1                 0.41                      0.57                0.043                 0.15 To numerically approximate y(1), where (1+t^2)y''+y'-y=0, y(0)=1,y'(0)=-1, using 4 steps of Euler's method, first convert to a system: y1' = y2, y1(0)=1; y2' = (y1-y2)/(1+t^2), y2(0)=-1. To numerically approximate y(1), where (1+t^2)y''+y'-y=0, y(0)=1, y'(0)=-1, using 4 steps of Euler's method, first convert to a system: y_1' = y_2, y_1(0)=1; y_2' = \\frac{y_1-y_2}{1+t^2}, y_2(0)=-1.:: sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: t, x, y=PolynomialRing(RR,3,"txy").gens() 3/4                 0.63                   -0.0078               -0.031                 0.11 1                 0.63                     0.020                0.079                0.071 To numerically approximate y(1), where y''+ty'+y=0, y(0)=1,y'(0)=0: To numerically approximate y(1), where y''+ty'+y=0, y(0)=1, y'(0)=0:: sage: t,x,y=PolynomialRing(RR,3,"txy").gens() sage: f = y; g = -x-y*t 3/4                 0.88                     -0.15                -0.62                -0.10 1                 0.75                     -0.17                -0.68               -0.015 AUTHOR: David Joyner AUTHORS: - David Joyner """ if method=="table": print "%10s %20s %25s %20s %20s"%("t", "x","h*f(t,x,y)","y", "h*g(t,x,y)") def eulers_method_2x2_plot(f,g, t0, x0, y0, h, t1): """ This plots the soln in the rectangle (xrange[0],xrange[1]) x (yrange[0],yrange[1]) and plots using Euler's method the numerical solution of the 1st order ODEs x' = f(t,x,y), x(a)=x0, y' = g(t,x,y), y(a) = y0. This plots the soln in the rectangle (xrange[0],xrange[1]) x (yrange[0],yrange[1]) and plots using Euler's method the numerical solution of the 1st order ODEs x' = f(t,x,y), x(a)=x_0, y' = g(t,x,y), y(a) = y_0. *For pedagogical purposes only.* === The following example plots the solution to theta''+sin(theta)=0, theta(0)=3/4, theta'(0) = 0.  Type P[0].show() to plot the solution, (P[0]+P[1]).show() to plot (t,theta(t)) and (t,theta'(t)). === EXAMPLES: EXAMPLES:: The following example plots the solution to \\theta''+\\sin(\\theta)=0, \\theta(0)=\\frac 34, \\theta'(0) = 0.  Type P[0].show() to plot the solution, (P[0]+P[1]).show() to plot (t,\\theta(t)) and (t,\\theta'(t)):: sage: from sage.calculus.desolvers import eulers_method_2x2_plot sage: f = lambda z : z[2]; g = lambda z : -sin(z[1]) sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0) """ Used to determine bounds for numerical integration. If end_points is None, the interval for integration is from ics[0] to ics[0]+10 - If end_points is None, the interval for integration is from ics[0] to ics[0]+10 If end_points is a or [a], the interval for integration is from min(ics[0],a) to max(ics[0],a) - If end_points is a or [a], the interval for integration is from min(ics[0],a) to max(ics[0],a) If end_points is [a,b], the interval for integration is from min(ics[0],a) to max(ics[0],b) - If end_points is [a,b], the interval for integration is from min(ics[0],a) to max(ics[0],b) EXAMPLES: EXAMPLES:: sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds sage: desolve_rk4_determine_bounds([0,2],1) (0, 1) :: sage: desolve_rk4_determine_bounds([0,2]) (0, 10) :: sage: desolve_rk4_determine_bounds([0,2],[-2]) (-2, 0) :: sage: desolve_rk4_determine_bounds([0,2],[-2,4]) (-2, 4) def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds): """ Solves numerically one first-order ordinary differential equation. Solves numerically one first-order ordinary differential equation. See also ode_solver. INPUT: input is similar to desolve command. The differential equation can be written in a form close to the plot_slope_field or desolve command input is similar to desolve command. The differential equation can be written in a form close to the plot_slope_field or desolve command Variant 1 (function in two variables) de   -- right hand side, i.e. the function f(x,y) from ODE y'=f(x,y) dvar -- dependent variable (symbolic variable declared by var) - Variant 1 (function in two variables) - de - right hand side, i.e. the function f(x,y) from ODE y'=f(x,y) - dvar - dependent variable (symbolic variable declared by var) Variant 2 (symbolic equation) de   -- equation, including term with diff(y,x) dvar -- dependent variable (declared as funciton of independent variable) - Variant 2 (symbolic equation) Other parameters ivar  -- should be specified, if there are more variables or if the equation is autonomous ics   -- initial conditions in the form [x0,y0] end_points  -- the end points of the interval if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) if end_points is None, we use end_points=ics[0]+10 if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) step   -- the length of the step (positive number) output -- 'list', 'plot', 'slope_field' (graph of the solution with slope field) - de - equation, including term with diff(y,x) - dvar - dependent variable (declared as funciton of independent variable) OUTPUT: Returns a list of points, or plot produced by list_plot, optionally with slope field. - Other parameters - ivar - should be specified, if there are more variables or if the equation is autonomous EXAMPLES: - ics - initial conditions in the form [x0,y0] - end_points - the end points of the interval - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is None, we use end_points=ics[0]+10 - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - step - (optional, default:0.1) the length of the step (positive number) - output - (optional, default: 'list') one of 'list', 'plot', 'slope_field' (graph of the solution with slope field) OUTPUT: Returns a list of points, or plot produced by list_plot, optionally with slope field. EXAMPLES:: sage: from sage.calculus.desolvers import desolve_rk4 Variant 2 for input - more common in numerics Variant 2 for input - more common in numerics:: sage: x,y=var('x y') sage: desolve_rk4(x*y*(2-y),y,ics=[0,1],end_points=1,step=0.5) Variant 1 for input - we can pass ODE in the form used by desolve function In this example we integrate bakwards, since end_points < ics[0] end_points < ics[0]:: sage: y=function('y',x) sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0) [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]] Here we show how to plot simple pictures. For more advanced aplications use list_plot instead. To see the resulting picture use show(P) in Sage notebook. Here we show how to plot simple pictures. For more advanced aplications use list_plot instead. To see the resulting picture use show(P) in Sage notebook. :: sage: x,y=var('x y') sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3) ALGORITHM: 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package. Perhaps could be faster by using fast_float instead. AUTHOR: Robert Marik (10-2009) 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.  Perhaps could be faster by using fast_float instead. AUTHORS: - Robert Marik (10-2009) """ if ics is None: raise ValueError, "No initial conditions, specify with ics=[x0,y0]." """ Solves numerically system of first-order ordinary differential equations using the 4th order Runge-Kutta method. Wrapper for Maxima command rk. Maxima command rk. See also ode_solver. INPUT: input is similar to desolve_system and desolve_rk4 commands des  -- right hand sides of the system vars -- dependent variables ivar -- should be specified, if there are more variables or if the equation is autonomous and the independent variable is missing ics  -- initial conditions in the form [x0,y01,y02,y03,....] end_points  -- the end points of the interval if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) if end_points is None, we use end_points=ics[0]+10 if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) step -- the length of the step input is similar to desolve_system and desolve_rk4 commands - des - right hand sides of the system - vars - dependent variables - ivar - (optional) should be specified, if there are more variables or if the equation is autonomous and the independent variable is missing - ics - initial conditions in the form [x0,y01,y02,y03,....] - end_points - the end points of the interval - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is None, we use end_points=ics[0]+10 - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - step -- (optional, default: 0.1) the length of the step OUTPUT: Returns a list of points. Returns a list of points. EXAMPLES: Lotka Volterra system EXAMPLES:: Lotka Volterra system:: sage: from sage.calculus.desolvers import desolve_system_rk4 sage: x,y,t=var('x y t') sage: LP=list_plot(Q) ALGORITHM: 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package. Perhaps could be faster by using fast_float instead. AUTHOR: Robert Marik (10-2009) 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.  Perhaps could be faster by using fast_float instead. AUTHOR: - Robert Marik (10-2009) """ if ics is None:
• ## sage/gsl/ode.pyx

diff -r b89cf68ccc1c -r ea16778cc7ae sage/gsl/ode.pyx
 a include 'gsl.pxi' import sage.gsl.interpolation import sage.gsl.interpolation cdef class PyFunctionWrapper: cdef object the_function cdef object the_jacobian cdef object the_parameters cdef int y_n cdef object the_function cdef object the_jacobian cdef object the_parameters cdef int y_n cdef set_yn(self,x): self.y_n = x cdef set_yn(self,x): self.y_n = x cdef class ode_system: cdef int  c_j(self,double t, double *y, double *dfdy,double *dfdt): #void *params): return 0 cdef int  c_j(self,double t, double *y, double *dfdy,double *dfdt): #void *params): return 0 cdef int c_f(self,double t, double* y, double* dydt):  #void *params): return 0 cdef int c_f(self,double t, double* y, double* dydt):  #void *params): return 0 cdef int c_jac_compiled(double t, double *y, double *dfdy,double *dfdt, void * params): cdef int status cdef ode_system wrapper wrapper = params status = wrapper.c_j(t,y,dfdy,dfdt)    #Could add parameters return status cdef int status cdef ode_system wrapper wrapper = params status = wrapper.c_j(t,y,dfdy,dfdt)    #Could add parameters return status cdef int c_f_compiled(double t, double *y, double *dydt, void *params): cdef int status cdef ode_system wrapper wrapper = params status =  wrapper.c_f(t,y,dydt)  #Could add parameters return status cdef int status cdef ode_system wrapper wrapper = params status =  wrapper.c_f(t,y,dydt)  #Could add parameters return status cdef int c_jac(double t,double *y,double *dfdy,double *dfdt,void *params): cdef int i cdef int j cdef int y_n cdef int param_n cdef PyFunctionWrapper wrapper wrapper = params y_n=wrapper.y_n y_list=[] for i from 0<=i params y_n=wrapper.y_n y_list=[] for i from 0<=i params y_n= wrapper.y_n y_list=[] for i from 0<=i params y_n= wrapper.y_n y_list=[] for i from 0<=i2: raise ValueError, "ODE system has a parameter but no parameters specified" elif self.params!=[]: wrapper.the_parameters = self.params wrapper.y_n = dim def plot_solution(self, i=0, filename=None, interpolate=False): from sage.plot.all import plot, point points=[] for x in self.solution: points.append(point((x[0],x[1][i]))) t = plot(points) if filename is None: t.show() else: t.save(filename=filename) def ode_solve(self,t_span=False,y_0=False,num_points=False,params=[]): import inspect cdef double h # step size h=self.h cdef int i cdef int j cdef int type cdef int dim cdef PyFunctionWrapper wrapper #struct to pass information into GSL C function self.params=params cdef double t cdef double t_end cdef double *y cdef double * scale_abs_array scale_abs_array=NULL y= malloc(sizeof(double)*(dim)) if y==NULL: raise MemoryError,"error allocating memory" result=[] v=[0]*dim cdef gsl_odeiv_step_type * T if t_span != False: self.t_span = t_span if y_0 != False: self.y_0 = y_0 for i from 0 <=i< dim: #copy initial conditions into C array y[i]=self.y_0[i] dim = len(self.y_0) type = isinstance(self.function,ode_system) if type == 0: wrapper = PyFunctionWrapper() if self.function!=None: wrapper.the_function = self.function else: raise ValueError, "ODE system not yet defined" if self.jacobian is None: wrapper.the_jacobian = None else: wrapper.the_jacobian = self.jacobian if self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])==2: wrapper.the_parameters=[] elif self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])>2: raise ValueError, "ODE system has a parameter but no parameters specified" elif self.params!=[]: wrapper.the_parameters = self.params wrapper.y_n = dim if self.algorithm == "rkf45": T=gsl_odeiv_step_rkf45 elif self.algorithm == "rk2": T=gsl_odeiv_step_rk2 elif self.algorithm == "rk4": T=gsl_odeiv_step_rk4 elif self.algorithm == "rkck": T=gsl_odeiv_step_rkck elif self.algorithm == "rk8pd": T=gsl_odeiv_step_rk8pd elif self.algorithm == "rk2imp": T= gsl_odeiv_step_rk2imp elif self.algorithm == "rk4imp": T= gsl_odeiv_step_rk4imp elif self.algorithm == "bsimp": T = gsl_odeiv_step_bsimp if not type and self.jacobian==None: raise TypeError,"The jacobian must be provided for the implicit Burlisch-Stoer method" elif self.algorithm == "gear1": T = gsl_odeiv_step_gear1 elif self.algorithm == "gear2": T = gsl_odeiv_step_gear2 else: raise TypeError,"algorithm not valid" cdef double t cdef double t_end cdef double *y cdef double * scale_abs_array scale_abs_array=NULL cdef gsl_odeiv_step * s s  = gsl_odeiv_step_alloc (T, dim) if s==NULL: free(y) raise MemoryError, "error setting up solver" y= malloc(sizeof(double)*(dim)) if y==NULL: raise MemoryError,"error allocating memory" result=[] v=[0]*dim cdef gsl_odeiv_step_type * T for i from 0 <=i< dim: #copy initial conditions into C array y[i]=self.y_0[i] cdef gsl_odeiv_control * c if self.algorithm == "rkf45": T=gsl_odeiv_step_rkf45 elif self.algorithm == "rk2": T=gsl_odeiv_step_rk2 elif self.algorithm == "rk4": T=gsl_odeiv_step_rk4 elif self.algorithm == "rkck": T=gsl_odeiv_step_rkck elif self.algorithm == "rk8pd": T=gsl_odeiv_step_rk8pd elif self.algorithm == "rk2imp": T= gsl_odeiv_step_rk2imp elif self.algorithm == "rk4imp": T= gsl_odeiv_step_rk4imp elif self.algorithm == "bsimp": T = gsl_odeiv_step_bsimp if not type and self.jacobian==None: raise TypeError,"The jacobian must be provided for the implicit Burlisch-Stoer method" elif self.algorithm == "gear1": T = gsl_odeiv_step_gear1 elif self.algorithm == "gear2": T = gsl_odeiv_step_gear2 else: raise TypeError,"algorithm not valid" if self.a == False and self.a_dydt==False: c  = gsl_odeiv_control_y_new (self.error_abs, self.error_rel) elif self.a !=False and self.a_dydt != False: if self.scale_abs==False: c = gsl_odeiv_control_standard_new(self.error_abs,self.error_rel,self.a,self.a_dydt) elif hasattr(self.scale_abs,'__len__'): if len(self.scale_abs)==dim: scale_abs_array = malloc(dim*sizeof(double)) for i from 0 <=i malloc(dim*sizeof(double)) for i from 0 <=iself.function).the_parameters = self.params sys.params = self.function else:                  # The user passed a python function. sys.function = c_f sys.jacobian = c_jac sys.params = wrapper sys.dimension = dim cdef int status import copy cdef int n if len(self.t_span)==2 and num_points!=False: try: n = num_points except TypeError: gsl_odeiv_evolve_free (e) if c == NULL: gsl_odeiv_control_free (c) gsl_odeiv_step_free (s) free(y) if scale_abs_array!=NULL: free(scale_abs_array) raise TypeError,"numpoints must be integer" result.append( (self.t_span[0],self.y_0)) delta = (self.t_span[1]-self.t_span[0])/(1.0*num_points) t =self.t_span[0] t_end=self.t_span[0]+delta for i from 0 y[j] result.append( (t,copy.copy(v)) ) t = t_end t_end= t+delta cdef gsl_odeiv_evolve * e e  = gsl_odeiv_evolve_alloc(dim) elif len(self.t_span) > 2: n = len(self.t_span) result.append((self.t_span[0],self.y_0)) t=self.t_span[0] t_end=self.t_span[1] for i from 0self.function).the_parameters = self.params sys.params = self.function else:                  # The user passed a python function. sys.function = c_f sys.jacobian = c_jac sys.params = wrapper sys.dimension = dim for j from 0<=j y[j] result.append( (t,copy.copy(v)) ) t=t_span[i] t_end=t_span[i+1] cdef int status import copy cdef int n if len(self.t_span)==2 and num_points!=False: try: n = num_points except TypeError: gsl_odeiv_evolve_free (e) gsl_odeiv_control_free (c) gsl_odeiv_step_free (s) free(y) if scale_abs_array!=NULL: free(scale_abs_array) raise TypeError,"numpoints must be integer" result.append( (self.t_span[0],self.y_0)) delta = (self.t_span[1]-self.t_span[0])/(1.0*num_points) t =self.t_span[0] t_end=self.t_span[0]+delta for i from 0 y[j] result.append( (t,copy.copy(v)) ) t = t_end t_end= t+delta elif len(self.t_span) > 2: n = len(self.t_span) result.append((self.t_span[0],self.y_0)) t=self.t_span[0] t_end=self.t_span[1] for i from 0 y[j] result.append( (t,copy.copy(v)) ) t=t_span[i] t_end=t_span[i+1] gsl_odeiv_evolve_free (e) gsl_odeiv_control_free (c) gsl_odeiv_step_free (s) free(y) if scale_abs_array!=NULL: free(scale_abs_array) self.solution = result
• ## sage/misc/functional.py

diff -r b89cf68ccc1c -r ea16778cc7ae sage/misc/functional.py
 a def integral(x, *args, **kwds): """ Returns an indefinite integral of an object x. Returns an indefinite or definite integral of an object x. First call x.integrate() and if that fails make an object and integrate it using Maxima, maple, etc, as specified by algorithm. For symbolic expression calls sage.calculus.calculus.integral` - see this function for available options. EXAMPLES:: sage: f = cyclotomic_polynomial(10) sage: integral(f) 1/5*x^5 - 1/4*x^4 + 1/3*x^3 - 1/2*x^2 + x :: sage: integral(sin(x),x) -cos(x) :: sage: y = var('y') sage: integral(sin(x),y) y*sin(x) :: sage: integral(sin(x), x, 0, pi/2) 1 sage: sin(x).integral(x, 0,pi/2)