Symbolic Math Toolbox

Solving Algebraic and Differential Equations

This example shows how to solve algebraic and differential equations using Symbolic Math Toolbox.

Solving Algebraic Equations

The solve function can handle a wide variety of algebraic equations, including both single equations and systems of equations. By default, the solver always tries to return a complete set of solutions, including complex solutions and multiple branches. For example, the following polynomial has two real and two complex solutions.

solve(x^4 + x^3 - x = 1, x)

{-1, 1, - 1/2 + (3^(1/2)*I)/2, - 1/2 - (3^(1/2)*I)/2}

The Real option lets you omit complex solutions and return only real ones:

solve(x^4 + x^3 - x = 1, x, Real)

{-1, 1}

By default, the solve function returns a very complete solution that takes into account all edge case situations.In situations where you know these edge cases will not occur, you can specify the IgnoreAnalyticConstraints option.This option switches the solver to a less strict mode of computation, and is very useful if you want to get a simpler solution that is typically valid but potentially incomplete for unlikely edge case situations.You can solve the same equation with and without this option to see the difference.

solve(ln(y) + ln(x + 1) = 1, y)

piecewise([Im(ln(x + 1)) in Dom::Interval([-PI], PI), {exp(1 - ln(x + 1))}], [not Im(ln(x + 1)) in Dom::Interval([-PI], PI), {}])

solve(ln(y) + ln(x + 1) = 1, y, IgnoreAnalyticConstraints)

{exp(1)/(x + 1)}

Some equations, such as those that include trigonometric elements, have an infinite number of solutions. For such equations, the solver returns the solution in the form of membership.

solve(sin(x) = 1, x)

Dom::ImageSet(PI/2 + 2*PI*k, k, Z_)

To get only one solution from the set, you can use the PrincipalValue option.

solve(sin(x) = 1, x, PrincipalValue)

{PI/2}

Solving Systems of Algebraic Equations

The solve function can also be used to solve systems of equations. Simply specify the system of equations and the variables to solve for.

solve([y + x^2 = 1, x - y = 10], [x, y], VectorFormat)

{matrix([[- (3*5^(1/2))/2 - 1/2], [- (3*5^(1/2))/2 - 21/2]]), matrix([[(3*5^(1/2))/2 - 1/2], [(3*5^(1/2))/2 - 21/2]])}

As with an individual equation, if a system has an infinite number of solutions, the solver returns the solution in the form of membership.

solve([exp(x + y) = 1, x - y = 1], [x, y], VectorFormat)

Dom::ImageSet(matrix([[1/2 + PI*k*I], [- 1/2 + PI*k*I]]), k, Z_)

In some cases, you may be able to express the system of equations in matrix form, A*`x→` = `b→`.For these cases, you can use the linalg::matlinsolve function to solve the equations.

A := matrix(2, 2, [1, 2, 3, 4]):
b := matrix(2, 1, [5, 6]):
x := linalg::matlinsolve(A, b)

matrix([[-4], [9/2]])

Solving Ordinary Differential Equations (ODEs)

You can solve various types of ordinary differential equations using the standard solve function (i.e. the same function that is used to solve algebraic equations) or the ode::solve function. If you use the standard solve function, you must first define your equation as an ODE.The same underlying algorithm is used regardless of which approach you choose.

You can solve a simple 2nd order linear ODE using the standard solve function.

eqn := ode(y''(t) + y'(t) + y(t) = 0, y(t))

ode((D@@2)(y)(t) + D(y)(t) + y(t), y(t))

solve(eqn)

{(C2*cos((3^(1/2)*t)/2))/exp(t/2) + (C3*sin((3^(1/2)*t)/2))/exp(t/2)}

You can solve the same equation using the ode::solve function, which does not require you to declare the equation as an ODE.

ode::solve(y''(t) + y'(t) + y(t) = 0, y(t))

{(C5*cos((3^(1/2)*t)/2))/exp(t/2) + (C6*sin((3^(1/2)*t)/2))/exp(t/2)}

You can specify initial conditions and boundary conditions along with an ODE, which is often critical for real world problems in engineering and physics.

f := ode::solve({y''(t) = 2*t/y'(t), y(0) = 0, y'(0) = 1}, y(t))

{(2^(1/2)*ln(t + (t^2 + 1/2)^(1/2)))/4 - (2^(1/2)*ln(2^(1/2)/2))/4 + (2^(1/2)*t*(t^2 + 1/2)^(1/2))/2}

You can visualize the solution to your differential equation using the plot function.Because the solver returns the result as a set, you must extract the solution from the set using brackets.

plot(f[1])

MuPAD graphics