Coverage for local/lib/python2.7/site-packages/sage/calculus/desolvers.py : 72%

r""" 

Solving ordinary differential equations 

This file contains functions useful for solving differential equations 

which occur commonly in a 1st semester differential equations 

course. For another numerical solver see the :meth:`ode_solver` function 

and the optional package Octave. 

Solutions from the Maxima package can contain the three constants 

``_C``, ``_K1``, and ``_K2`` where the underscore is used to distinguish 

them from symbolic variables that the user might have used. You can 

substitute values for them, and make them into accessible usable 

symbolic variables, for example with ``var("_C")``. 

Commands: 

- :func:`desolve` - Compute the "general solution" to a 1st or 2nd order 

ODE via Maxima. 

- :func:`desolve_laplace` - Solve an ODE using Laplace transforms via 

Maxima. Initial conditions are optional. 

- :func:`desolve_rk4` - Solve numerically an IVP for one first order 

equation, return list of points or plot. 

- :func:`desolve_system_rk4` - Solve numerically an IVP for a system of first 

order equations, return list of points. 

- :func:`desolve_odeint` - Solve numerically a system of first-order ordinary 

differential equations using ``odeint`` from `scipy.integrate module. 

<https://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate>`_ 

- :func:`desolve_system` - Solve a system of 1st order ODEs of any size using 

Maxima. Initial conditions are optional. 

- :func:`eulers_method` - Approximate solution to a 1st order DE, 

presented as a table. 

- :func:`eulers_method_2x2` - Approximate solution to a 1st order system 

of DEs, presented as a table. 

- :func:`eulers_method_2x2_plot` - Plot the sequence of points obtained 

from Euler's method. 

The following functions require the optional package ``tides``: 

- :func:`desolve_mintides` - Numerical solution of a system of 1st order ODEs via 

the Taylor series integrator method implemented in TIDES. 

- :func:`desolve_tides_mpfr` - Arbitrary precision Taylor series integrator implemented in TIDES. 

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 

- Miguel Marco (06-2014) - Tides desolvers 

""" 

########################################################################## 

# Copyright (C) 2006 David Joyner <wdjoyner@gmail.com>, Marshall Hampton, 

# Robert Marik <marik@mendelu.cz> 

# 

# Distributed under the terms of the GNU General Public License (GPL): 

# 

# http://www.gnu.org/licenses/ 

########################################################################## 

from __future__ import division 

from sage.interfaces.maxima import Maxima 

from sage.plot.all import line 

from sage.symbolic.expression import is_SymbolicEquation 

from sage.symbolic.ring import is_SymbolicVariable 

from sage.calculus.functional import diff 

from sage.misc.functional import N 

from sage.misc.decorators import rename_keyword 

from tempfile import mkdtemp 

import shutil 

import os 

from sage.rings.real_mpfr import RealField 

maxima = Maxima() 

def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False): 

r""" 

Solves a 1st or 2nd order linear ODE via Maxima, including IVP and BVP. 

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`, i.e. write `[x_0, y(x_0), y'(x_0)]` 

- for a second-order boundary solution, specify initial and 

final `x` and `y` boundary conditions, i.e. write `[x_0, y(x_0), x_1, y(x_1)]`. 

- gives an error if the solution is not SymbolicEquation (as happens for 

example for a Clairaut 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 

the method which has been used to get a 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 the example of an equation which is separable but 

this property is not recognized by Maxima and the equation is solved 

as exact. 

- ``contrib_ode`` - (optional) if true, ``desolve`` allows to solve 

Clairaut, Lagrange, Riccati and some other equations. This may take 

a long time and is thus turned off by default. Initial conditions 

can be used only if the result is one SymbolicEquation (does not 

contain a singular solution, for example). 

OUTPUT: 

In most cases return a SymbolicEquation which defines the solution 

implicitly. If the result is in the form `y(x)=\ldots` (happens for 

linear eqs.), return the right-hand side only. The possible 

constant solutions of separable ODEs are omitted. 

NOTES: 

Use ``desolve? <tab>`` if the output in the Sage notebook is truncated. 

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) 

Graphics object consisting of 1 graphics primitive 

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) 

_K2*e^(-x) + _K1*e^x - x 

:: 

sage: f = desolve(de, y, [10,2,1]); f 

-x + 7*e^(x - 10) + 5*e^(-x + 10) 

:: 

sage: f(x=10) 

2 

:: 

sage: diff(f,x)(x=10) 

1 

:: 

sage: de = diff(y,x,2) + y == 0 

sage: desolve(de, y) 

_K2*cos(x) + _K1*sin(x) 

:: 

sage: desolve(de, y, [0,1,pi/2,4]) 

cos(x) + 4*sin(x) 

:: 

sage: desolve(y*diff(y,x)+sin(x)==0,y) 

-1/2*y(x)^2 == _C - cos(x) 

Clairaut 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 an 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]] 

:: 

sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True,show_method=True) 

[[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati'] 

Higher order equations, 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]).expand() 

1/6*y(x)^3 - 5/3*y(x) == x - 3/2 

:: 

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:: 

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,show_method=True) 

[-cos(y(x)) == _C + sin(x), 'separable'] 

:: 

sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1]) 

-cos(y(x)) == -cos(1) + sin(x) - 1 

Linear equation - Sage returns the expression on the right hand side only:: 

sage: desolve(diff(y,x)+(y) == cos(x),y) 

1/2*((cos(x) + sin(x))*e^x + 2*_C)*e^(-x) 

:: 

sage: desolve(diff(y,x)+(y) == cos(x),y,show_method=True) 

[1/2*((cos(x) + sin(x))*e^x + 2*_C)*e^(-x), 'linear'] 

:: 

sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1]) 

1/2*(cos(x)*e^x + e^x*sin(x) + 1)*e^(-x) 

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, also using initial 

conditions, but you cannot put (sometimes desired) the initial 

condition at `x=0`, since this point is a singular point of the 

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) 

Example of difficult ODE producing an 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." 

Another difficult ODE with error - moreover, it takes a long time :: 

sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested 

Some more types of ODEs:: 

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 an error (as expected, Maxima 5.18 cannot 

solve equations from initial conditions). 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() # 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:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y) 

(_K2*x + _K1)*e^(-x) + 1/2*sin(x) 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True) 

[(_K2*x + _K1)*e^(-x) + 1/2*sin(x), 'variationofparameters'] 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1]) 

1/2*(7*x + 6)*e^(-x) + 1/2*sin(x) 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True) 

[1/2*(7*x + 6)*e^(-x) + 1/2*sin(x), 'variationofparameters'] 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2]) 

3*(x*(e^(1/2*pi) - 2)/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*(x*(e^(1/2*pi) - 2)/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) 

[(_K2*x + _K1)*e^(-x), 'constcoeff'] 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1]) 

(4*x + 3)*e^(-x) 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True) 

[(4*x + 3)*e^(-x), 'constcoeff'] 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2]) 

(2*x*(2*e^(1/2*pi) - 3)/pi + 3)*e^(-x) 

:: 

sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True) 

[(2*x*(2*e^(1/2*pi) - 3)/pi + 3)*e^(-x), 'constcoeff'] 

TESTS: 

:trac:`9961` fixed (allow assumptions on the dependent variable in desolve):: 

sage: y=function('y')(x); assume(x>0); assume(y>0) 

sage: sage.calculus.calculus.maxima('domain:real') # needed since Maxima 5.26.0 to get the answer as below 

real 

sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True) 

[x - arcsinh(y(x)/x) == _C] 

:trac:`10682` updated Maxima to 5.26, and it started to show a different 

solution in the complex domain for the ODE above:: 

sage: forget() 

sage: sage.calculus.calculus.maxima('domain:complex') # back to the default complex domain 

complex 

sage: assume(x>0) 

sage: assume(y>0) 

sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True) 

[x - arcsinh(y(x)^2/(x*sqrt(y(x)^2))) - arcsinh(y(x)/x) + 1/2*log(4*(x^2 + 2*y(x)^2 + 2*sqrt(x^2*y(x)^2 + y(x)^4))/x^2) == _C] 

: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 

:trac:`9835` fixed:: 

sage: x = var('x') 

sage: y = function('y')(x) 

sage: desolve(diff(y,x,2)+y*(1-y^2)==0,y,[0,-1,1,1]) 

Traceback (most recent call last): 

... 

NotImplementedError: Unable to use initial condition for this equation (freeofx). 

:trac:`8931` fixed:: 

sage: x=var('x'); f=function('f')(x); k=var('k'); assume(k>0) 

sage: desolve(diff(f,x,2)/f==k,f,ivar=x) 

_K1*e^(sqrt(k)*x) + _K2*e^(-sqrt(k)*x) 

:trac:`15775` fixed:: 

sage: forget() 

sage: y = function('y')(x) 

sage: desolve(diff(y, x) == sqrt(abs(y)), dvar=y, ivar=x) 

sqrt(-y(x))*(sgn(y(x)) - 1) + (sgn(y(x)) + 1)*sqrt(y(x)) == _C + x 

AUTHORS: 

- David Joyner (1-2006) 

- Robert Bradshaw (10-2008) 

- Robert Marik (10-2009) 

""" 

if is_SymbolicEquation(de): 

de = de.lhs() - de.rhs() 

if is_SymbolicVariable(dvar): 

raise ValueError("You have to declare dependent variable as a function evaluated at the independent variable, eg. y=function('y')(x)") 

# for backwards compatibility 

if isinstance(dvar, list): 

dvar, ivar = dvar 

elif ivar is None: 

ivars = de.variables() 

ivars = [t for t in ivars if t is not dvar] 

if len(ivars) != 1: 

raise ValueError("Unable to determine independent variable, please specify.") 

ivar = ivars[0] 

de00 = de._maxima_() 

P = de00.parent() 

dvar_str=P(dvar.operator()).str() 

ivar_str=P(ivar).str() 

de00 = de00.str() 

def sanitize_var(exprs): 

return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) 

de0 = sanitize_var(de00) 

ode_solver="ode2" 

cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str) 

# we produce string like this 

# ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x) 

soln = P(cmd) 

if str(soln).strip() == 'false': 

if contrib_ode: 

ode_solver="contrib_ode" 

P("load('contrib_ode)") 

cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str) 

# we produce string like this 

# (TEMP:contrib_ode(x*('diff(y,x,1))^2-(x*y+1)*'diff(y,x,1)+y,y,x), if TEMP=false then TEMP else substitute(y=y(x),TEMP)) 

soln = P(cmd) 

if str(soln).strip() == 'false': 

raise NotImplementedError("Maxima was unable to solve this ODE.") 

else: 

raise NotImplementedError("Maxima was unable to solve this ODE. Consider to set option contrib_ode to True.") 

if show_method: 

maxima_method=P("method") 

if (ics is not None): 

if not is_SymbolicEquation(soln.sage()): 

if not show_method: 

maxima_method=P("method") 

raise NotImplementedError("Unable to use initial condition for this equation (%s)."%(str(maxima_method).strip())) 

if len(ics) == 2: 

tempic=(ivar==ics[0])._maxima_().str() 

tempic=tempic+","+(dvar==ics[1])._maxima_().str() 

cmd="(TEMP:ic1(%s(%s,%s,%s),%s),substitute(%s=%s(%s),TEMP))"%(ode_solver,de00,dvar_str,ivar_str,tempic,dvar_str,dvar_str,ivar_str) 

cmd=sanitize_var(cmd) 

# we produce string like this 

# (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP)) 

soln=P(cmd) 

if len(ics) == 3: 

#fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 

P("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \ 

noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \ 

temp: lhs(soln) - rhs(soln), \ 

TEMP_k:solve([subst([xa,ya],soln), subst([dya,xa], lhs(dya)=-subst(0,lhs(dya),diff(temp,lhs(xa)))/diff(temp,lhs(ya)))],[%k1,%k2]), \ 

if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \ 

temp: maplist(lambda([zz], subst(zz,soln)), TEMP_k), \ 

if length(temp)=1 then return(first(temp)) else return(temp))") 

tempic=P(ivar==ics[0]).str() 

tempic=tempic+","+P(dvar==ics[1]).str() 

tempic=tempic+",'diff("+dvar_str+","+ivar_str+")="+P(ics[2]).str() 

cmd="(TEMP:ic2_sage(%s(%s,%s,%s),%s),substitute(%s=%s(%s),TEMP))"%(ode_solver,de00,dvar_str,ivar_str,tempic,dvar_str,dvar_str,ivar_str) 

cmd=sanitize_var(cmd) 

# we produce string like this 

# (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP)) 

soln=P(cmd) 

if str(soln).strip() == 'false': 

raise NotImplementedError("Maxima was unable to solve this IVP. Remove the initial condition to get the general solution.") 

if len(ics) == 4: 

#fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 

P("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \ 

noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \ 

TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \ 

if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \ 

temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \ 

if length(temp)=1 then return(first(temp)) else return(temp))") 

cmd="bc2_sage(%s(%s,%s,%s),%s,%s=%s,%s,%s=%s)"%(ode_solver,de00,dvar_str,ivar_str,P(ivar==ics[0]).str(),dvar_str,P(ics[1]).str(),P(ivar==ics[2]).str(),dvar_str,P(ics[3]).str()) 

cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str) 

cmd=sanitize_var(cmd) 

# we produce string like this 

# (TEMP:bc2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,x=%pi/2,y=2),substitute(y=y(x),TEMP)) 

soln=P(cmd) 

if str(soln).strip() == 'false': 

raise NotImplementedError("Maxima was unable to solve this BVP. Remove the initial condition to get the general solution.") 

soln=soln.sage() 

if is_SymbolicEquation(soln) and soln.lhs() == dvar: 

# Remark: Here we do not check that the right hand side does not depend on dvar. 

# This probably will not happen for solutions obtained via ode2, anyway. 

soln = soln.rhs() 

if show_method: 

return [soln,maxima_method.str()] 

else: 

return soln 

#def desolve_laplace2(de,vars,ics=None): 

## """ 

## Solves an ODE using laplace transforms via maxima. Initial conditions 

## are optional. 

## INPUT: 

## de -- a lambda expression representing the ODE 

## (eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)") 

## vars -- a list of strings representing the variables 

## (eg, vars = ["x","f"], if x is the independent 

## variable and f is the dependent variable) 

## ics -- a list of numbers representing initial conditions, 

## with symbols allowed which are represented by strings 

## (eg, f(0)=1, f'(0)=2 is ics = [0,1,2]) 

## EXAMPLES:: 

## sage: from sage.calculus.desolvers import desolve_laplace 

## sage: x = var('x') 

## sage: f = function('f')(x) 

## sage: de = lambda y: diff(y,x,x) - 2*diff(y,x) + y 

## sage: desolve_laplace(de(f(x)),[f,x]) 

## #x*%e^x*(?%at('diff('f(x),x,1),x=0))-'f(0)*x*%e^x+'f(0)*%e^x 

## sage: desolve_laplace(de(f(x)),[f,x],[0,1,2]) ## IC option does not work 

## #x*%e^x*(?%at('diff('f(x),x,1),x=0))-'f(0)*x*%e^x+'f(0)*%e^x 

## AUTHOR: David Joyner (1st version 1-2006, 8-2007) 

## """ 

# ######## this method seems reasonable but doesn't work for some reason 

# name0 = vars[0]._repr_()[0:(len(vars[0]._repr_())-2-len(str(vars[1])))] 

# name1 = str(vars[1]) 

# #maxima("de:"+de+";") 

# if ics is not None: 

# ic0 = maxima("ic:"+str(vars[1])+"="+str(ics[0])) 

# d = len(ics) 

# for i in range(d-1): 

# maxima(vars[0](vars[1])).diff(vars[1],i).atvalue(ic0,ics[i+1]) 

# de0 = de._maxima_() 

# #cmd = "desolve("+de+","+vars[1]+"("+vars[0]+"));" 

# #return maxima.eval(cmd) 

# return de0.desolve(vars[0]).rhs() 

def desolve_laplace(de, dvar, ics=None, ivar=None): 

""" 

Solve an ODE using Laplace transforms. Initial conditions are optional. 

INPUT: 

- ``de`` - a lambda expression representing the ODE (e.g. ``de = 

diff(y,x,2) == diff(y,x)+sin(x)``) 

- ``dvar`` - the dependent variable (e.g. ``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, (e.g. 

``f(0)=1``, ``f'(0)=2`` corresponds to ``ics = [0,1,2]``) 

OUTPUT: 

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:: 

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):: 

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 

sage: desolve_laplace(de,f) 

-x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0) 

:: 

sage: desolve_laplace(de,f,ics=[0,1,2]) 

x*e^x + e^x 

TESTS: 

Check that :trac:`4839` is fixed:: 

sage: t = var('t') 

sage: x = function('x')(t) 

sage: soln = desolve_laplace(diff(x,t)+x==1, x, ics=[0,2]) 

sage: soln 

e^(-t) + 1 

:: 

sage: soln(t=3) 

e^(-3) + 1 

AUTHORS: 

- 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;") 

#if ics is not None: 

# d = len(ics) 

# for i in range(0,d-1): 

# ic = "atvalue(diff("+vars[1]+"("+vars[0]+"),"+str(vars[0])+","+str(i)+"),"+str(vars[0])+"="+str(ics[0])+","+str(ics[1+i])+")" 

# maxima(ic) 

# 

#cmd = "desolve("+de._repr_()+","+vars[1]+"("+vars[0]+"));" 

#return maxima(cmd).rhs()._maxima_init_() 

## verbatim copy from desolve - begin 

if is_SymbolicEquation(de): 

de = de.lhs() - de.rhs() 

if is_SymbolicVariable(dvar): 

raise ValueError("You have to declare dependent variable as a function evaluated at the independent variable, eg. y=function('y')(x)") 

# for backwards compatibility 

if isinstance(dvar, list): 

dvar, ivar = dvar 

elif ivar is None: 

ivars = de.variables() 

ivars = [t for t in ivars if t is not dvar] 

if len(ivars) != 1: 

raise ValueError("Unable to determine independent variable, please specify.") 

ivar = ivars[0] 

## verbatim copy from desolve - end 

dvar_str = str(dvar) 

def sanitize_var(exprs): # 'y(x) -> y(x) 

return exprs.replace("'"+dvar_str,dvar_str) 

de0=de._maxima_() 

P = de0.parent() 

i = dvar_str.find('(') 

dvar_str = dvar_str[:i+1] + '_SAGE_VAR_' + dvar_str[i+1:] 

cmd = sanitize_var("desolve("+de0.str()+","+dvar_str+")") 

soln=P(cmd).rhs() 

if str(soln).strip() == 'false': 

raise NotImplementedError("Maxima was unable to solve this ODE.") 

soln=soln.sage() 

if ics is not None: 

d = len(ics) 

for i in range(0,d-1): 

soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') 

return soln 

def desolve_system(des, vars, ics=None, ivar=None): 

""" 

Solve a system of any size of 1st order ODEs. Initial conditions are optional. 

One dimensional systems are passed to :meth:`desolve_laplace`. 

INPUT: 

- ``des`` - list of ODEs 

- ``vars`` - list of dependent variables 

- ``ics`` - (optional) list of initial values for ``ivar`` and ``vars``. 

If ``ics`` is defined, it should provide initial conditions for each variable, 

otherwise an exception would be raised. 

- ``ivar`` - (optional) the independent variable, which must be 

specified if there is more than one independent variable in the 

equation. 

EXAMPLES:: 

sage: t = var('t') 

sage: x = function('x')(t) 

sage: y = function('y')(t) 

sage: de1 = diff(x,t) + y - 1 == 0 

sage: de2 = diff(y,t) - x + 1 == 0 

sage: desolve_system([de1, de2], [x,y]) 

[x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1, 

y(t) == (y(0) - 1)*cos(t) + (x(0) - 1)*sin(t) + 1] 

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)) # not tested 

sage: parametric_plot((solnx,solny),(0,1)) # not tested 

TESTS: 

Check that :trac:`9823` is fixed:: 

sage: t = var('t') 

sage: x = function('x')(t) 

sage: de1 = diff(x,t) + 1 == 0 

sage: desolve_system([de1], [x]) 

-t + x(0) 

Check that :trac:`16568` is fixed:: 

sage: t = var('t') 

sage: x = function('x')(t) 

sage: y = function('y')(t) 

sage: de1 = diff(x,t) + y - 1 == 0 

sage: de2 = diff(y,t) - x + 1 == 0 

sage: des = [de1,de2] 

sage: ics = [0,1,-1] 

sage: vars = [x,y] 

sage: sol = desolve_system(des, vars, ics); sol 

[x(t) == 2*sin(t) + 1, y(t) == -2*cos(t) + 1] 

:: 

sage: solx, soly = sol[0].rhs(), sol[1].rhs() 

sage: RR(solx(t=3)) 

1.28224001611973 

:: 

sage: P1 = plot([solx,soly], (0,1)) 

sage: P2 = parametric_plot((solx,soly), (0,1)) 

Now type ``show(P1)``, ``show(P2)`` to view these plots. 

Check that :trac:`9824` is fixed:: 

sage: t = var('t') 

sage: epsilon = var('epsilon') 

sage: x1 = function('x1')(t) 

sage: x2 = function('x2')(t) 

sage: de1 = diff(x1,t) == epsilon 

sage: de2 = diff(x2,t) == -2 

sage: desolve_system([de1, de2], [x1, x2], ivar=t) 

[x1(t) == epsilon*t + x1(0), x2(t) == -2*t + x2(0)] 

sage: desolve_system([de1, de2], [x1, x2], ics=[1,1], ivar=t) 

Traceback (most recent call last): 

... 

ValueError: Initial conditions aren't complete: number of vars is different from number of dependent variables. Got ics = [1, 1], vars = [x1(t), x2(t)] 

AUTHORS: 

- Robert Bradshaw (10-2008) 

- Sergey Bykov (10-2014) 

""" 

if ics is not None: 

if len(ics) != (len(vars) + 1): 

raise ValueError("Initial conditions aren't complete: number of vars is different from number of dependent variables. Got ics = {0}, vars = {1}".format(ics, vars)) 

if len(des)==1: 

return desolve_laplace(des[0], vars[0], ics=ics, ivar=ivar) 

ivars = set([]) 

for i, de in enumerate(des): 

if not is_SymbolicEquation(de): 

des[i] = de == 0 

ivars = ivars.union(set(de.variables())) 

if ivar is None: 

ivars = ivars - set(vars) 

if len(ivars) != 1: 

raise ValueError("Unable to determine independent variable, please specify.") 

ivar = list(ivars)[0] 

dvars = [v._maxima_() for v in vars] 

if ics is not None: 

ivar_ic = ics[0] 

for dvar, ic in zip(dvars, ics[1:]): 

dvar.atvalue(ivar==ivar_ic, ic) 

soln = dvars[0].parent().desolve(des, dvars) 

if str(soln).strip() == 'false': 

raise NotImplementedError("Maxima was unable to solve this system.") 

soln = list(soln) 

for i, sol in enumerate(soln): 

soln[i] = sol.sage() 

if ics is not None: 

ivar_ic = ics[0] 

for dvar, ic in zip(dvars, ics[:1]): 

dvar.atvalue(ivar==ivar_ic, dvar) 

return soln 

@rename_keyword(deprecation=6094, method="algorithm") 

def eulers_method(f,x0,y0,h,x1,algorithm="table"): 

r""" 

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 `x_0` to `x_1` by `h` (so 

`(x_1-x_0)/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. 

.. NOTE:: 

This function is for pedagogical purposes only. 

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) 

x y h*f(x,y) 

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,algorithm="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,algorithm="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) 

x y h*f(x,y) 

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) 

1 1 1/3 

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,algorithm="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,algorithm="none") 

sage: P1 = list_plot(pts) 

sage: P2 = line(pts) 

sage: (P1+P2).show() 

AUTHORS: 

- David Joyner 

""" 

if algorithm=="table": 

print("%10s %20s %25s"%("x","y","h*f(x,y)")) 

n=int((1.0)*(x1-x0)/h) 

x00=x0; y00=y0 

soln = [[x00,y00]] 

for i in range(n+1): 

if algorithm=="table": 

print("%10r %20r %20r"%(x00,y00,h*f(x00,y00))) 

y00 = y00+h*f(x00,y00) 

x00=x00+h 

soln.append([x00,y00]) 

if algorithm!="table": 

return soln 

@rename_keyword(deprecation=6094, method="algorithm") 

def eulers_method_2x2(f,g, t0, x0, y0, h, t1,algorithm="table"): 

r""" 

This implements Euler's method for finding numerically the 

solution of the 1st order system of two ODEs 

.. MATH:: 

\begin{aligned} 

x' &= f(t, x, y), x(t_0)=x_0 \\ 

y' &= g(t, x, y), y(t_0)=y_0. 

\end{aligned} 

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. 

.. NOTE:: 

This function is for pedagogical purposes only. 

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,algorithm="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 

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.00 0 0.00 

1/3 0.00 0.13 0.00 0.00 

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: `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() 

sage: f = y; g = (x-y)/(1+t^2) 

sage: eulers_method_2x2(f,g, 0, 1, -1, 1/4, 1) 

t x h*f(t,x,y) y h*g(t,x,y) 

0 1 -0.25 -1 0.50 

1/4 0.75 -0.12 -0.50 0.29 

1/2 0.63 -0.054 -0.21 0.19 

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`:: 

sage: t,x,y=PolynomialRing(RR,3,"txy").gens() 

sage: f = y; g = -x-y*t 

sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1) 

t x h*f(t,x,y) y h*g(t,x,y) 

0 1 0.00 0 -0.25 

1/4 1.0 -0.062 -0.25 -0.23 

1/2 0.94 -0.11 -0.46 -0.17 

3/4 0.88 -0.15 -0.62 -0.10 

1 0.75 -0.17 -0.68 -0.015 

AUTHORS: 

- David Joyner 

""" 

if algorithm=="table": 

print("%10s %20s %25s %20s %20s"%("t", "x","h*f(t,x,y)","y", "h*g(t,x,y)")) 

n=int((1.0)*(t1-t0)/h) 

t00 = t0; x00 = x0; y00 = y0 

soln = [[t00,x00,y00]] 

for i in range(n+1): 

if algorithm=="table": 

print("%10r %20r %25r %20r %20r"%(t00,x00,h*f(t00,x00,y00),y00,h*g(t00,x00,y00))) 

x01 = x00 + h*f(t00,x00,y00) 

y00 = y00 + h*g(t00,x00,y00) 

x00 = x01 

t00 = t00 + h 

soln.append([t00,x00,y00]) 

if algorithm!="table": 

return soln 

def eulers_method_2x2_plot(f,g, t0, x0, y0, h, t1): 

r""" 

Plot solution of ODE. 

This plots the solution in the rectangle with sides ``(xrange[0],xrange[1])`` and 

``(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`. 

.. NOTE:: 

This function is for pedagogical purposes only. 

EXAMPLES: