Solving Systems of ODEs and Linear Systems with Sage

[Eli]

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

• 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).

Output

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

Code: [Select all] [Expand/Collapse]

1. x = var('x')
2. y = function('y')(x)
3. desolve(diff(y,x) + y - 1, y)
4.  
5. #Plot
6. x = var('x')
7. y = function('y')(x)
8. f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f
9. plot(f)

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(y*diff(y,x)+sin(x)==0,y)

Second order differential equations

Code: [Select all] [Expand/Collapse]

1. x = var('x')
2. y = function('y')(x)
3. de = diff(y,x,2) - y == x
4. desolve(de, y)
5.  
6. f = desolve(de, y, [10,2,1]); f
7. plot(f)

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. de = diff(y,x,2) + y == 0
3. desolve(de, y)

Code: [Select all] [Expand/Collapse]

1. #BVP
2. y = function('y')(x)
3. de = diff(y,x,2) + y == 0
4. desolve(de, y, [0,1,pi/2,4])

Higher order equations, not involving independent variable

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
3.  
4. #Solution 2
5. y = function('y')(x)
6. desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()
7.  
8. #Solution 3
9. y = function('y')(x)
10. 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

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(diff(y,x)*sin(y) == cos(x),y)
3.  
4. #Solution 2
5. y = function('y')(x)
6. desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True)
7.  
8. #Solution 3
9. y = function('y')(x)
10. desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])

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

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(diff(y,x)+(y) == cos(x),y)
3.  
4. #Solution 2
5. y = function('y')(x)
6. desolve(diff(y,x)+(y) == cos(x),y,show_method=True)
7.  
8. #Solution 3
9. y = function('y')(x)
10. f = desolve(diff(y,x)+(y) == cos(x),y,[0,1]); f
11. 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}\):

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. 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.

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. f = desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y); f

Some more types of ODEs

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
3.  
4. #Solution 2
5. y = function('y')(x)
6. desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)

Second order linear ODE

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
3.  
4. #Solution 2
5. y = function('y')(x)
6. desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True)
7.  
8. #Solution 3
9. y = function('y')(x)
10. desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1])
11.  
12. #Solution 4
13. y = function('y')(x)
14. 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:

Code: [Select all] [Expand/Collapse]

1. y = function('y')(x)
2. de = x^3*diff(y, x, 3) + x^2*diff(y, x, 2) - 2*x*diff(y, x) + 2*y - 2*x^4
3. Y = desolve(de, y, algorithm="fricas"); Y  #Optional - fricas
4.  
5. #Solution 2
6.  
7. y = function('y')(x)
8. de = x^3*diff(y, x, 3) + x^2*diff(y, x, 2) - 2*x*diff(y, x) + 2*y - 2*x^4
9. 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)

Code: [Select all] [Expand/Collapse]

1. t = var('t')
2. x = function('x')(t)
3. y = function('y')(t)
4. de1 = diff(x,t) + y - 1 == 0
5. de2 = diff(y,t) - x + 1 == 0
6. desolve_system([de1, de2], [x,y])
7.  
8. #Solve system of ODEs using FriCAS
9. desolve_system([de1, de2], [x,y], algorithm='fricas')
10.  
11. #Apply some initial conditions
12. 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:

Code: [Select all] [Expand/Collapse]

1. from sage.calculus.desolvers import desolve_odeint
2. x,y=var('x,y')
3. f=[x*(1-y),-y*(1-x)]
4. sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])
5. p=line(zip(sol[:,0],sol[:,1]))
6. p.show()

Lorenz Equations:

Code: [Select all] [Expand/Collapse]

1. x,y,z=var('x,y,z')
2. # Next we define the parameters
3. sigma=10
4. rho=28
5. beta=8/3
6. # The Lorenz equations
7. lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
8. # Time and initial conditions
9. times=srange(0,50.05,0.05)
10. ics=[0,1,1]
11. 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:

Code: [Select all] [Expand/Collapse]

1. from sage.calculus.desolvers import desolve_system_rk4
2. x,y,t=var('x y t')
3. P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
4. Q=[ [i,j] for i,j,k in P]
5. LP=list_plot(Q)
6.  
7. Q=[ [j,k] for i,j,k in P]
8. LP=list_plot(Q)