## Solve System of Algebraic Equations

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.

### Handle the Output of solve

Suppose you have the system

`$\begin{array}{c}{x}^{2}{y}^{2}=0\\ x-\frac{y}{2}=\alpha ,\end{array}$`

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`

### Solve a Linear System of Equations

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.

### Return the Full Solution of a System of Equations

`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:

`$\begin{array}{l}\mathrm{sin}\left(x\right)\text{​}+\mathrm{cos}\left(y\right)=\frac{4}{5}\\ \mathrm{sin}\left(x\right)\mathrm{cos}\left(y\right)=\frac{1}{10}\end{array}$`

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'))```

### Solve a System of Equations Under Conditions

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')``` ### Work with Solutions, Parameters, and Conditions Returned by solve

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``` ### Convert Symbolic Results to Numeric Values

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]```

### Simplify Complicated Results and Improve Performance

If results look complicated, `solve` is stuck, or if you want to improve performance, see, Troubleshoot Equation Solutions from solve Function.

## Support

#### Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos