Solving ordinary differential equations
See references at SageMath Website.
desolve
sage.calculus.desolvers.desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False, algorithm='maxima')
Computes the general solution of the first or second order linear ODEs, including IVP and BVP, via Maxima.
Inputs
- 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 [x0,y(x0),y'(x0)],
- for a second-order boundary solution, specify initial and final x and y boundary conditions, i.e. write [x0,y(x0),x1,y(x1)],
- 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).
- algorithm – (default: 'maxima') one of:
- 'maxima' - use maxima,
- 'fricas' - use FriCAS (the optional fricas spkg has to be installed).
In most cases return a SymbolicEquation which defines the solution implicitly. If the result is in the form y(x)=… (happens for linear eqs.), return the right-hand side only. The possible constant solutions of separable ODEs are omitted.
desolve usage
First order differential equations
- x = var('x')
- y = function('y')(x)
- desolve(diff(y,x) + y - 1, y)
- #Plot
- x = var('x')
- y = function('y')(x)
- f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f
- plot(f)
- y = function('y')(x)
- desolve(y*diff(y,x)+sin(x)==0,y)
Second order differential equations
- x = var('x')
- y = function('y')(x)
- de = diff(y,x,2) - y == x
- desolve(de, y)
- f = desolve(de, y, [10,2,1]); f
- plot(f)
- y = function('y')(x)
- de = diff(y,x,2) + y == 0
- desolve(de, y)
- #BVP
- y = function('y')(x)
- de = diff(y,x,2) + y == 0
- desolve(de, y, [0,1,pi/2,4])
Higher order equations, not involving independent variable
- y = function('y')(x)
- desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
- #Solution 2
- y = function('y')(x)
- desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()
- #Solution 3
- y = function('y')(x)
- desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
Separable equations - Sage returns solution in implicit form
- y = function('y')(x)
- desolve(diff(y,x)*sin(y) == cos(x),y)
- #Solution 2
- y = function('y')(x)
- desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True)
- #Solution 3
- y = function('y')(x)
- desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
Linear equation - Sage returns the expression on the right hand side only
- y = function('y')(x)
- desolve(diff(y,x)+(y) == cos(x),y)
- #Solution 2
- y = function('y')(x)
- desolve(diff(y,x)+(y) == cos(x),y,show_method=True)
- #Solution 3
- y = function('y')(x)
- f = desolve(diff(y,x)+(y) == cos(x),y,[0,1]); f
- plot(f)
This ODE with separated variables is solved as exact. Explanation - factor does not split \(e^{x−y}\) in Maxima into \(e^xe^{-y}\):
- y = function('y')(x)
- desolve(diff(y,x)==exp(x-y),y,show_method=True)
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. If the solution should be bounded at \(x=0\), then _K_2=0.
- y = function('y')(x)
- f = desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y); f
Some more types of ODEs
- y = function('y')(x)
- desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
- #Solution 2
- y = function('y')(x)
- desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
Second order linear ODE
- y = function('y')(x)
- desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
- #Solution 2
- y = function('y')(x)
- desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True)
- #Solution 3
- y = function('y')(x)
- desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1])
- #Solution 4
- y = function('y')(x)
- desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True)
Using algorithm='fricas' we can invoke the differential equation solver from FriCAS. For example, it can solve higher order linear equations:
- y = function('y')(x)
- de = x^3*diff(y, x, 3) + x^2*diff(y, x, 2) - 2*x*diff(y, x) + 2*y - 2*x^4
- Y = desolve(de, y, algorithm="fricas"); Y #Optional - fricas
- #Solution 2
- y = function('y')(x)
- de = x^3*diff(y, x, 3) + x^2*diff(y, x, 2) - 2*x*diff(y, x) + 2*y - 2*x^4
- Y = desolve(de, y, ics=[1,3,7], algorithm="fricas"); Y #Optional - fricas
Here, the initial conditions are then interpreted as \([(x0, \ y(x0), \ y'(x0),…,y^{(n)}(x0)]\)
Solving Systems of ODEs
sage.calculus.desolvers.desolve_system(des, vars, ics=None, ivar=None, algorithm='maxima')
Solves a system of any size of 1st order ODEs. Initial conditions are optional.
One dimensional systems are passed to 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.
algorithm – (default: 'maxima') one of:
- 'maxima' - use maxima
- 'fricas' - use FriCAS (the optional fricas spkg has to be installed)
- t = var('t')
- x = function('x')(t)
- y = function('y')(t)
- de1 = diff(x,t) + y - 1 == 0
- de2 = diff(y,t) - x + 1 == 0
- desolve_system([de1, de2], [x,y])
- #Solve system of ODEs using FriCAS
- desolve_system([de1, de2], [x,y], algorithm='fricas')
- #Apply some initial conditions
- sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol
sage.calculus.desolvers.desolve_odeint(des, ics, times, dvars, ivar=None, compute_jac=False, args=(), rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
Solves numerically a system of first-order ordinary differential equations using odeint from scipy.integrate module.
Return a list with the solution of the system at each time in times.
Lotka Volterra Equations:
- from sage.calculus.desolvers import desolve_odeint
- x,y=var('x,y')
- f=[x*(1-y),-y*(1-x)]
- sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])
- p=line(zip(sol[:,0],sol[:,1]))
- p.show()
Lorenz Equations:
- x,y,z=var('x,y,z')
- # Next we define the parameters
- sigma=10
- rho=28
- beta=8/3
- # The Lorenz equations
- lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
- # Time and initial conditions
- times=srange(0,50.05,0.05)
- ics=[0,1,1]
- sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)
sage.calculus.desolvers.desolve_system_rk4(des, vars, ics=None, ivar=None, end_points=None, step=0.1)
Solves numerically a system of first-order ordinary differential equations using the 4th order Runge-Kutta method. Wrapper for Maxima command rk.
Return a list of points.
Lotka Volterra system:
- from sage.calculus.desolvers import desolve_system_rk4
- x,y,t=var('x y t')
- P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
- Q=[ [i,j] for i,j,k in P]
- LP=list_plot(Q)
- Q=[ [j,k] for i,j,k in P]
- LP=list_plot(Q)