x′y′==−5xx+−3y3y(ODE)
x(0)=0 and y(0)=4 (ics).
# Phase Plane var('x y'); l = 5 xp = -5*x + 3*y yp = x - 3*y xnull = implicit_plot(xp == 0, (x,-l,l), (y,-l,l), color='red') ynull = implicit_plot(yp == 0, (x,-l,l), (y,-l,l), color='blue') # The first four initial conditions yield the straight line solutions. # The other initial conditions were chosen to illustrate all solution types. ics = [(1,1), (-1,-1), (3,-1), (-3,1), (0,4), (0,-4), (3,0), (-3,0)] solns = streamline_plot([xp, yp], (x,-l,l), (y,-l,l), start_points = ics, color='green') pics = point(ics) show(xnull + ynull + solns + pics, aspect_ratio=1)
# Set up the matrix of coefficients and vector of variables A = matrix([[-5, 3], [1, -3]]) var('x y'); X = vector([x, y]) show(A, X, A*X)
# Equilibria eqns = [0 == (A*X)[k] for k in range(2)]; show(eqns) solve(eqns, [x,y]) # A \ vector([0, 0]) returns (0, 0) regardless of A and so is not helpful.
# Alternatively, we know that the zero vector is an equilibrium. # This will be the unique equilibrium if the determinent of A is nonzero. A.det()
# With the guess X = B*e^(r*t), we obtain AB = rB, or equivalently (A - rI)B = 0. # This will have a non-zero solution B iff det(A - rI) = 0. # det(A - rI) is called the characteristic polynomial of A. var('r') AmrI = A - r*identity_matrix(2); show(AmrI) cpol = AmrI.det(); show(cpol); show(cpol.expand()) solve(cpol, r)
# The solutions to det(A - rI) = 0 are called the eigenvalues of A. A.eigenvalues()
# For r = -6, we solve (A - rI)B = 0 for the B vector to obtain (SOLN1). var('a b') RHS = (A - (-6)*identity_matrix(2))*vector([a,b]) eqns = [0 == RHS[k] for k in range(2)]; show(eqns) solve(eqns, [a,b])
# For r = -2, we solve (A - rI)B = 0 for the B vector to obtain (SOLN2). var('a b') RHS = (A - (-2)*identity_matrix(2))*vector([a,b]) eqns = [0 == RHS[k] for k in range(2)]; show(eqns) solve(eqns, [a,b])
# The eigenvalues and vectors can be obtained in one step. (evalues, evectors) = A.eigenmatrix_left(); show(evalues, evectors)
# The general solution can be created from the eigen information. var('a b t') xsol(t) = a*1*e^(-2*t) + b* 3 *e^(-6*t); show(xsol) ysol(t) = a*1*e^(-2*t) + b*(-1)*e^(-6*t); show(ysol)
# Find the coefficients for the specific solution. solve([xsol(0) == 0, ysol(0) == 4], [a,b])
# Find the general solution. var('t') x = function('x')(t) y = function('y')(t) de1 = diff(x,t) == -5*x + 3*y de2 = diff(y,t) == x - 3*y sol = desolve_system([de1, de2], [x, y], ivar = t); print(sol); show(sol)
# Use the general solution to find the specific solution. xsol(t) = sage_eval(str(sol[0].rhs()).replace('x(0)', '0').replace('y(0)', '4'), locals = {'t':t}); show(xsol) ysol(t) = sage_eval(str(sol[1].rhs()).replace('x(0)', '0').replace('y(0)', '4'), locals = {'t':t}); show(ysol)
# Find the specific solution directly. sol = desolve_system([de1, de2], [x, y], ivar = t, ics = [0, 0, 4]); show(sol) xsol(t) = sol[0].rhs(); show(xsol) ysol(t) = sol[1].rhs(); show(ysol)
Note that it is sometimes useful to go backward, as well as forward, in time.
# Can you tell which color is x(t) and which color is y(t)? plot([xsol, ysol], (-0.1,2), figsize = 4, axes_labels = ['t','x,y'])
# Where is (x(0), y(0))? (x(2), y(2))? parametric_plot([xsol, ysol], (-0.1, 2))