This topic shows you how to solve a system of equations symbolically using Symbolic Math Toolbox™. This toolbox offers both numeric and symbolic equation solvers. For a comparison of numeric and symbolic solvers, see Select Numeric or Symbolic Solver.
Suppose you have the system
and you want to solve for x and y. First, create the necessary symbolic objects.
syms x y a
There are several ways to address the output of solve
.
One way is to use a two-output call.
[solx,soly] = solve(x^2*y^2 == 0, x-y/2 == a)
The call returns the following.
solx = 0 a soly = -2*a 0
Modify the first equation to x2y2 = 1. The new system has more solutions.
[solx,soly] = solve(x^2*y^2 == 1, x-y/2 == a)
Four distinct solutions are produced.
solx = a/2 - (a^2 - 2)^(1/2)/2 a/2 - (a^2 + 2)^(1/2)/2 a/2 + (a^2 - 2)^(1/2)/2 a/2 + (a^2 + 2)^(1/2)/2 soly = - a - (a^2 - 2)^(1/2) - a - (a^2 + 2)^(1/2) (a^2 - 2)^(1/2) - a (a^2 + 2)^(1/2) - a
Since you did not specify the dependent variables, solve
uses symvar
to
determine the variables.
This way of assigning output from solve
is
quite successful for “small” systems. For instance,
if you have a 10-by-10 system of equations, typing the following is
both awkward and time consuming.
[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10] = solve(...)
To circumvent this difficulty, solve
can
return a structure whose fields are the solutions. For example, solve
the system of equations u^2 - v^2 = a^2
, u
+ v = 1
, a^2 - 2*a = 3
.
syms u v a S = solve(u^2 - v^2 == a^2, u + v == 1, a^2 - 2*a == 3)
The solver returns its results enclosed in this structure.
S = struct with fields: a: [2×1 sym] u: [2×1 sym] v: [2×1 sym]
The solutions for a
reside in the “a
-field”
of S
.
S.a
ans = -1 3
Similar comments apply to the solutions for u
and v
.
The structure S
can now be manipulated by the field
and index to access a particular portion of the solution. For example,
to examine the second solution, you can use the following statement
to extract the second component of each field.
s2 = [S.a(2), S.u(2), S.v(2)]
s2 = [ 3, 5, -4]
The following statement creates the solution matrix M
whose
rows comprise the distinct solutions of the system.
M = [S.a, S.u, S.v]
M = [ -1, 1, 0] [ 3, 5, -4]
Clear solx
and soly
for
further use.
clear solx soly
Linear systems of equations can also be solved using matrix division. For example, solve this system.
clear u v x y syms u v x y eqns = [x + 2*y == u, 4*x + 5*y == v]; S = solve(eqns); sol = [S.x; S.y] [A,b] = equationsToMatrix(eqns,x,y); z = A\b
sol = (2*v)/3 - (5*u)/3 (4*u)/3 - v/3 z = (2*v)/3 - (5*u)/3 (4*u)/3 - v/3
Thus,sol
and z
produce
the same solution, although the results are assigned to different
variables.
solve
does not automatically return all
solutions of an equation. To return all solutions along with the parameters
in the solution and the conditions on the solution, set the ReturnConditions
option
to true
.
Consider the following system of equations:
Visualize the system of equations using fimplicit
.
To set the x-axis and y-axis
values in terms of pi
, get the axes handles using axes
in a
.
Create the symbolic array S
of the values -2*pi
to 2*pi
at
intervals of pi/2
. To set the ticks to S
,
use the XTick
and YTick
properties
of a
. To set the labels for the x-and y-axes, convert S
to
character vectors. Use arrayfun
to apply char
to
every element of S
to return T
.
Set the XTickLabel
and YTickLabel
properties
of a
to T
.
syms x y eqn1 = sin(x)+cos(y) == 4/5; eqn2 = sin(x)*cos(y) == 1/10; a = axes; fimplicit(eqn1,[-2*pi 2*pi],'b'); hold on grid on fimplicit(eqn2,[-2*pi 2*pi],'m'); L = sym(-2*pi:pi/2:2*pi); a.XTick = double(L); a.YTick = double(L); M = arrayfun(@char, L, 'UniformOutput', false); a.XTickLabel = M; a.YTickLabel = M; title('Plot of System of Equations') legend('sin(x)+cos(y) == 4/5','sin(x)*cos(y) == 1/10',... 'Location','best','AutoUpdate','off')
The solutions lie at the intersection of the two plots. This
shows the system has repeated, periodic solutions. To solve this system
of equations for the full solution set, use solve
and
set the ReturnConditions
option to true
.
S = solve(eqn1, eqn2, 'ReturnConditions', true)
S = struct with fields: x: [2×1 sym] y: [2×1 sym] parameters: [1×2 sym] conditions: [2×1 sym]
solve
returns a structure S
with
the fields S.x
for the solution to x
, S.y
for
the solution to y
, S.parameters
for
the parameters in the solution, and S.conditions
for
the conditions on the solution. Elements of the same index in S.x
, S.y
,
and S.conditions
form a solution. Thus, S.x(1)
, S.y(1)
,
and S.conditions(1)
form one solution to the system
of equations. The parameters in S.parameters
can
appear in all solutions.
Index into S
to return the solutions, parameters,
and conditions.
S.x S.y S.parameters S.conditions
ans = z1 z1 ans = z z ans = [ z, z1] ans = (in((z - acos(6^(1/2)/10 + 2/5))/(2*pi), 'integer') |... in((z + acos(6^(1/2)/10 + 2/5))/(2*pi), 'integer')) &... (in(-(pi - z1 + asin(6^(1/2)/10 - 2/5))/(2*pi), 'integer') |... in((z1 + asin(6^(1/2)/10 - 2/5))/(2*pi), 'integer')) (in((z1 - asin(6^(1/2)/10 + 2/5))/(2*pi), 'integer') |... in((z1 - pi + asin(6^(1/2)/10 + 2/5))/(2*pi), 'integer')) &... (in((z - acos(2/5 - 6^(1/2)/10))/(2*pi), 'integer') |... in((z + acos(2/5 - 6^(1/2)/10))/(2*pi), 'integer'))
To solve the system of equations under conditions, specify the
conditions in the input to solve
.
Solve the system of equations considered above for x
and y
in
the interval -2*pi
to 2*pi
.
Overlay the solutions on the plot using scatter
.
Srange = solve(eqn1, eqn2, -2*pi<x, x<2*pi, -2*pi<y, y<2*pi, 'ReturnConditions', true); scatter(Srange.x, Srange.y,'k')
You can use the solutions, parameters, and conditions returned
by solve
to find solutions within an interval
or under additional conditions. This section has the same goal as
the previous section, to solve the system of equations within a search
range, but with a different approach. Instead of placing conditions
directly, it shows how to work with the parameters and conditions
returned by solve
.
For the full solution S
of the system of
equations, find values of x
and y
in
the interval -2*pi
to 2*pi
by
solving the solutions S.x
and S.y
for
the parameters S.parameters
within that interval
under the condition S.conditions
.
Before solving for x
and y
in
the interval, assume the conditions in S.conditions
using assume
so
that the solutions returned satisfy the condition. Assume the conditions
for the first solution.
assume(S.conditions(1))
Find the parameters in S.x
and S.y
.
paramx = intersect(symvar(S.x), S.parameters) paramy = intersect(symvar(S.y), S.parameters)
paramx = z1 paramy = z
Solve the first solution of x
for the parameter paramx
.
solparamx(1,:) = solve(S.x(1) > -2*pi, S.x(1) < 2*pi, paramx)
solparamx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)]
Similarly, solve the first solution of y
for paramy
.
solparamy(1,:) = solve(S.y(1) > -2*pi, S.y(1) < 2*pi, paramy)
solparamy = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)]
Clear the assumptions set by S.conditions(1)
using assume
.
Call asumptions
to check that the assumptions
are cleared.
assume(S.parameters,'clear') assumptions
ans = Empty sym: 1-by-0
Assume the conditions for the second solution.
assume(S.conditions(2))
Solve the second solution to x
and y
for
the parameters paramx
and paramy
.
solparamx(2,:) = solve(S.x(2) > -2*pi, S.x(2) < 2*pi, paramx) solparamy(2,:) = solve(S.y(2) > -2*pi, S.y(2) < 2*pi, paramy)
solparamx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)] [ asin(6^(1/2)/10 + 2/5), pi - asin(6^(1/2)/10 + 2/5), asin(6^(1/2)/10 + 2/5) - 2*pi, - pi - asin(6^(1/2)/10 + 2/5)] solparamy = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)] [ acos(2/5 - 6^(1/2)/10), acos(2/5 - 6^(1/2)/10) - 2*pi, -acos(2/5 - 6^(1/2)/10), 2*pi - acos(2/5 - 6^(1/2)/10)]
The first rows of paramx
and paramy
form
the first solution to the system of equations, and the second rows
form the second solution.
To find the values of x
and y
for
these values of paramx
and paramy
,
use subs
to substitute for paramx
and paramy
in S.x
and S.y
.
solx(1,:) = subs(S.x(1), paramx, solparamx(1,:)); solx(2,:) = subs(S.x(2), paramx, solparamx(2,:)) soly(1,:) = subs(S.y(1), paramy, solparamy(1,:)); soly(2,:) = subs(S.y(2), paramy, solparamy(2,:))
solx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)] [ asin(6^(1/2)/10 + 2/5), pi - asin(6^(1/2)/10 + 2/5), asin(6^(1/2)/10 + 2/5) - 2*pi, - pi - asin(6^(1/2)/10 + 2/5)] soly = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)] [ acos(2/5 - 6^(1/2)/10), acos(2/5 - 6^(1/2)/10) - 2*pi, -acos(2/5 - 6^(1/2)/10), 2*pi - acos(2/5 - 6^(1/2)/10)]
Note that solx
and soly
are
the two sets of solutions to x
and to y
.
The full sets of solutions to the system of equations are the two
sets of points formed by all possible combinations of the values in solx
and soly
.
Plot these two sets of points using scatter
.
Overlay them on the plot of the equations. As expected, the solutions
appear at the intersection of the plots of the two equations.
for i = 1:length(solx(1,:)) for j = 1:length(soly(1,:)) scatter(solx(1,i), soly(1,j), 'k') scatter(solx(2,i), soly(2,j), 'k') end end
Symbolic calculations provide exact accuracy, while numeric
calculations are approximations. Despite this loss of accuracy, you
might need to convert symbolic results to numeric approximations for
use in numeric calculations. For a high-accuracy conversion, use variable-precision
arithmetic provided by the vpa
function.
For standard accuracy and better performance, convert to double precision
using double
.
Use vpa
to convert the symbolic solutions solx
and soly
to
numeric form.
vpa(solx) vpa(soly)
ans = [ 2.9859135500977407388300518406219,... -3.2972717570818457380952349259371,... 0.15567910349205249963259154265761,... -6.1275062036875339772926952239014] ... [ 0.70095651347102524787213653614929,... 2.4406361401187679905905068471302,... -5.5822287937085612290531502304097,... -3.8425491670608184863347799194288] ans = [ 0.86983981332387137135918515549046,... -5.4133454938557151055661016110685,... -0.86983981332387137135918515549046,... 5.4133454938557151055661016110685] ... [ 1.4151172233028441195987301489821,... -4.8680680838767423573265566175769,... -1.4151172233028441195987301489821,... 4.8680680838767423573265566175769]
If results look complicated, solve
is stuck,
or if you want to improve performance, see, Troubleshoot Equation Solutions from solve Function.