Solving Systems of ODEs and Linear Systems with Sage
Posted: Mon Jan 31, 2022 8:11 am
Run your computations here
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
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
Second order differential equations
Higher order equations, not involving independent variable
Separable equations - Sage returns solution in implicit form
Linear equation - Sage returns the expression on the right hand side only
This ODE with separated variables is solved as exact. Explanation - factor does not split \(e^{x−y}\) in Maxima into \(e^xe^{-y}\):
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.
Some more types of ODEs
Second order linear ODE
Using algorithm='fricas' we can invoke the differential equation solver from FriCAS. For example, it can solve higher order linear equations:
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:
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:
Lorenz Equations:
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:
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)