Today, you will learn about solving systems of linear equations using Sage.
{{{id=6| /// }}}You can create a system of linear equations that literally looks like a system of equations, then solve it for all the variables using the solve command.
{{{id=225| var('x,y') solve([2*x + 3*y == 5, 10*x + sqrt(2)*y == 1], [x,y]) /// [[x == -36/223*sqrt(2) + 35/446, y == 24/223*sqrt(2) + 360/223]] }}}NOTE (!): Above, we use == (double equals) to create symbolic equations. When you accidentally put = (single equals), you will get a syntax error.
In the following example, there are infinitely many solutions:
{{{id=213| var('x,y') solve([2*x + 3*y == 5, 4*x + 6*y == 10], [x,y]) /// [[x == -3/2*r2 + 5/2, y == r2]] }}}The main advantage of this approach is that it is so similar to the mathematics. However, it is cumbersome and potentially slow if the system is large.
Another advantage is that exactly the same format works for non-linear equations:
{{{id=230| var('x,y') solve([2*x^2 + 3*y^2 == 5, 4*x^2 + 6*y^3 == 10], [x,y]) /// [[x == -1/2*sqrt(2)*sqrt(5), y == 0], [x == 1/2*sqrt(2)*sqrt(5), y == 0], [x == -1, y == 1], [x == 1, y == 1]] }}}We plot the zero loci of the two equations above, along with the set of intersection points, as computed by the solve command. Note the use of the solution_dict=True parameter to obtain this solutions in an easier to use format.
{{{id=227| var('x,y') S = solve([2*x^2 + 3*y^2 == 5, 4*x^2 + 6*y^3 == 10], [x,y], solution_dict=True) print S G = implicit_plot(2*x^2 + 3*y^2 == 5, (x,-2,2), (y,-2,2)) G += implicit_plot(4*x^2 + 6*y^3 == 10, (x,-2,2), (y,-2,2), cmap='hsv') G += points([(a[x], a[y]) for a in S], pointsize=50) G /// [{y: 0, x: -1/2*sqrt(2)*sqrt(5)}, {y: 0, x: 1/2*sqrt(2)*sqrt(5)}, {y: 1, x: -1}, {y: 1, x: 1}]
}}}
{{{id=232|
///
}}}
{{{id=236|
///
}}}
Problem: Use the Sage solve command to solve the following system of three equations in three unknowns:
$2x + 3y + 5z = 7$
$11x +13 y +17 z = 19$
$23x + 29y + 31z = 37$
Problem: Use the Sage solve command to the following system of equations (note: Sage can't find exact solutions, so it returns numerical ones):
$2x^3 + \sqrt{2} \cdot y = x$
$x^2 + y^2 = 5$
In a linear algebra class you learn how to encode a system of linear equations as a matrix equation:
$$ A \cdot x = b $$
where $A$ is a matrix, $x$ is an unknown vector, and $b$ is a vector. Sage has excellent support for solving a system of linear equations presented this way.
Recall:
{{{id=244| var('x,y') solve([2*x + 3*y == 5, 10*x + sqrt(2)*y == 1], [x,y]) /// [[x == -36/223*sqrt(2) + 35/446, y == 24/223*sqrt(2) + 360/223]] }}}In terms of matrices:
{{{id=248| A = matrix(2, 2, [2,3, 10, sqrt(2)]); A /// [ 2 3] [ 10 sqrt(2)] }}} {{{id=243| b = vector([5, 1]); b /// (5, 1) }}} {{{id=250| x = A.solve_right(b); x /// (36/(sqrt(2) - 15) + 5/2, -24/(sqrt(2) - 15)) }}} {{{id=249| A*x == b /// True }}}Instead of solve_right, you can write A \ b like in MATLAB.
{{{id=242| A \ b /// (36/(sqrt(2) - 15) + 5/2, -24/(sqrt(2) - 15)) }}}Note that solve_right just finds one solution. It gives no information about how many solutions there are.
{{{id=253| var('x,y') solve([2*x + 3*y == 5, 4*x + 6*y == 10], [x,y]) /// [[x == -3/2*r3 + 5/2, y == r3]] }}} {{{id=255| A = matrix(QQ, 2, 2, [2,3, 4,6]) b = vector([5,10]) A.solve_right(b) /// (5/2, 0) }}}If you know linear algebra, you'll know that you can get all solutions by finding one solution, and adding to it any element of the right kernel of $A$. This right kernel is the set of vector $v$ such that $Av = 0$.
{{{id=252| A.right_kernel() /// Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 -2/3] }}} {{{id=257| var('r') * vector([1,-2/3]) + vector([5/2, 0]) /// (r + 5/2, -2/3*r) }}} {{{id=258| /// }}}Efficiency: If the matrix A has entries in the rational numbers or double precision floating point (e.g., RDF or CDF), then solving $Ax = b$ should be (not "is") very fast. In fact, Sage is one of the fastest programs in the world for doing this. Let's try some benchmarks.
A 50x50 system gets solved almost instantly:
{{{id=262| n = 50 A = random_matrix(ZZ, n) b = random_matrix(ZZ, n, 1) time x = A.solve_right(b) /// Time: CPU 0.02 s, Wall: 0.02 s }}}A 500x500 system only takes a few seconds:
{{{id=261| n = 500 A = random_matrix(ZZ, n) b = random_matrix(ZZ, n, 1) time x = A.solve_right(b) /// Time: CPU 2.61 s, Wall: 2.57 s }}}Writing down the solution vector takes about 1.5 million characters !
{{{id=263| len(x.str()) /// 1416999 }}}Here's the first coordinate, which is a huge rational number.
{{{id=259| x[0] /// (17149814220083302445442045813082967030641980151320029116122376006943203923985641215694936623436146037126669536219450456501249090830583284617332078645233312157181814080915962925465113013126437775068618601303857682717108861754827790031933857821646624177982513259124803316185656717300563399376865496128371187989193900657565416800593118564218666080298619474060325496458790156863813038023796282610711693830672126744648813978014311312732099382080001861358618170338386742431115051953116917620144396787088288614520038032102500036088960621742312904169206865981240712644415783840177379106728596940124598237111878191204932171135072067706120182654289498847508982357277720431056468504736560321189633550693453217159101725445305071751586711053045619629987468580529337922002043160668595993475226478027748488774841164570930836033916497682309125995636242750235148111579030177310450034900620546550383671932581877780276924718054389963784669541741268099661342647631027065594812437005151650256496771639868599741070213595701518955551123586530181750984030225142591138467897796520247358754680215067694748326999817706541866028421281726270087872423662850337267057443166564536794941973375616480022166983629145249016273102513535734397609297228674637217799823616913531102502263516392397532565230818851636490970852710627038151487610509308240321574221921953377014936544411649304660356086803485390010025896037531034643051044563447284474113170913199/9516176173938657193572938740758488962708232735912385179687953047925386627050062075805683377394408447459261662278880801127908978656435016266794546290874015800221145170001810367800941216621242595642006647601133340825458666292785270145232739846599701883875588679716465321719593542810492368045459063059666165855767491257961252222263997943301007808019609960478905389760595589807779944450638256112917344202994433132944767608476741681593961099579674888211504685612726171251725491491037930809181825907845787204420144152615403702611775639976303227406281140921509298132268468090613743835137432205472669208057340992734675119157992749053463309437899306014758502458199911118123779146289645050793515803648322720812353226316380738171515700111458477550178511161833192616748617831881106978145957218397064565680691626631203011215133214891152605790558776169292709767199180388603589002545552474425425242520743655700256130011886904648629998133736875400521753957605659625871171101595581094337220568344505938071006836018065356727026313208831231177123656690475736673470706008381023081032569131338914095058600606022207489288952448536675899400583765391388794341494105290886266739612065443855760669978449124733703732269470247354982372748022139015705317433473276562483672691280895669557223600124501901598541743974843041678034544915360670666892964039444011036045791127866291121008763241718765164222447544800863911551678541850005768901701352950) }}}Why so fast?
Here is a numerical example:
{{{id=270| n = 50 A = random_matrix(RDF, n) b = random_matrix(RDF, n, 1) time x = A.solve_right(b) /// Time: CPU 0.21 s, Wall: 0.22 s }}} {{{id=216| n = 200 A = random_matrix(RDF, n) b = random_matrix(RDF, n, 1) time x = A.solve_right(b) /// Time: CPU 11.95 s, Wall: 12.06 s }}}Yuck, that is DOG SLOW, because evidently nobody has implemented a fast solve_right using numpy yet (bummer!). We'll talk a lot more about numpy on Friday.
{{{id=215| import numpy time x = matrix(numpy.linalg.solve(A.numpy(), b.numpy())) /// Time: CPU 0.00 s, Wall: 0.00 s }}}Sanity check:
{{{id=212| (A*x - b)[:10] /// [-5.36376498772e-14] [-2.12052597703e-14] [-1.04916075827e-14] [ 1.13575815419e-13] [ 5.72875080707e-14] [-1.57651669497e-14] [ 2.93098878501e-14] [ 5.80091530367e-14] [ 2.16493489802e-14] [-8.65973959208e-15] }}} {{{id=221| n = 500 A = random_matrix(RDF, n) b = random_matrix(RDF, n, 1) import numpy time x = matrix(numpy.linalg.solve(A.numpy(), b.numpy())) /// Time: CPU 0.04 s, Wall: 0.03 s }}} {{{id=222| /// }}}Problem: Use Sage matrices to solve the following system of three equations in three unknowns:
$2x + 3y + 5z = 7$
$11x +13 y +17 z = 19$
$23x + 29y + 31z = 37$
A vector space is a a mathematical object -- like a ring -- that you will probably know about if you've taken linear algebra.
In Sage, every vector space is just the set of linear combinations of a finite list of vectors.
Sage fully supports vector spaces as mathematical objects in their own right.
Vector spaces arise natural as the kernels of matrices:
{{{id=224| A = random_matrix(QQ, 2, 4); A /// [-1/2 -1 2 -2] [ 1 -1 2 1] }}} {{{id=223| V = A.right_kernel(); V /// Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 -1/4 -1/2] [ 0 1 1/2 0] }}} {{{id=209| vector([0,2,1,0]) in V /// True }}} {{{id=278| V.basis() /// [ (1, 0, -1/4, -1/2), (0, 1, 1/2, 0) ] }}} {{{id=279| /// }}}