Documentation Center

  • Trial Software
  • Product Updates

Solve Equations Numerically

Get Numeric Results

There are two methods to get numeric approximations of the solutions:

  • Solve equations symbolically and approximate the obtained symbolic results numerically. Using this method, you get numeric approximations of all the solutions found by the symbolic solver. If the symbolic solver fails to find any solutions, MuPAD® calls the numeric solver directly. For nonpolynomial equations, the numeric solver returns only one solution. Using the symbolic solver and post-processing its results method requires more time than a purely numeric solver and can significantly decrease performance.

  • Solve equations using numeric methods from the beginning. This method is faster, but for nonpolynomial equations the numeric solver returns only the first solution it finds. You can use the AllRealRoots option to make the solver look for other real roots. Even with this option, the solution set can be incomplete. Using AllRealRoots can significantly decrease performance.

Approximate Symbolic Solutions Numerically

When you solve an equation symbolically, the solver does not always return results in an explicit form. For example, the solver can represent solutions using RootOf:

solve(x^3 + x + 1, x)

If you have a symbolic solution that you want to approximate numerically, use the float command:

float(%)

Use float to approximate the solutions of the symbolic system:

float(solve([x^3 + x^2 + 2*x = y, y^2 = x^2], [x, y]))

Approximating symbolic solutions numerically, you get the complete set of solutions. For example, solve the following equation symbolically:

S := solve(sin(x^2) = 1/2, x)

Suppose, you want to get numeric results instead of expressions containing PI. The float command returns the infinite solution set:

float(S)

Solve Equations Numerically

To avoid getting symbolic solutions and to proceed with numeric methods, use the numeric::solve command. If an equation is not polynomial, the numeric solver returns only one solution:

numeric::solve(sin(x^2) = 1/2, x)

Solve Polynomial Equations and Systems

For polynomial equations, the numeric solver returns all solutions:

numeric::solve(4*x^4 + 3*x^3 + 2*x^2 + x + 5 = 0, x)

When called with the option AllRealRoots, the solver omits all complex roots. For example, when you symbolically solve the polynomial equation and approximate the solutions, you get all solutions:

numeric::solve(4*x^4 + 3*x^3 + 2*x^2 + x -1 = 0, x)

To limit the solution set to the real solutions only, use the option AllRealRoots:

numeric::solve(4*x^4 + 3*x^3 + 2*x^2 + x - 1 = 0, x,
                                       AllRealRoots)

Using numeric::solve, you also can solve a system of polynomial equations:

numeric::solve([x^3 + 2*x = y, y^2 = x], [x, y])

To solve linear systems numerically, use the numeric::linsolve command. For example, solve the following system symbolically and numerically:

linsolve([x = y - 1, x + y = 5/2], [x, y]);
numeric::linsolve([x = y - 1, x + y = 5/2], [x, y])

Solve Arbitrary Algebraic Equations and Systems

For nonpolynomial equations, there is no general method of finding all the solutions numerically. When you solve a nonpolynomial equation or a system numerically, and the solutions exist, the solver returns only one solution:

numeric::solve(sin(1/x) = x, x)

This equation obviously has more than one solution:

plot(sin(1/x), x, x = -1..1)

To get more real solutions of a single equation containing one variable, call the numeric solver with the option AllRealRoots. The AllRealRoots option does not guarantee that the solver finds all existing real roots. For example, the option helps to find additional solutions for the equation:

numeric::solve(sin(1/x) = x, x, AllRealRoots)
Warning: Problem in isolating search intervals. Some roots might be lost. [numeric::allRealRoots]

For a system of nonpolynomial equation, the solver also returns only one solution. Plotting the equations, you see that the system has more than one solution:

numeric::solve([sin(x) = y^2 - 1, cos(x) = y], [x, y]);
plot(sin(x) = y^2 - 1, cos(x) = y)

The AllRealRoots option does not work for systems:

numeric::solve([sin(x) = y^2 - 1, cos(x) = y], [x, y],
                                         AllRealRoots)
Error: Only one equation is allowed with the 'AllRealRoots' option. [numeric::solve]

To find numerical approximations of other solutions, specify intervals that contain the solutions. You can use the command numeric::solve that internally calls numeric::fsolve. However, to speed up your calculations, call numeric::fsolve directly. Note that numeric::solve returns a set of solutions, and numeric::fsolve returns a list:

numeric::solve([sin(x) = y^2 - 1, cos(x) = y],
               [x = 2.5..3.5, y = -1.5..-0.5]);
numeric::fsolve([sin(x) = y^2 - 1, cos(x) = y],
                     [x = 4..5, y = -0.2..0.2])

The MultiSolutions option also serves to find more than one numeric approximation. Without this option, the numeric solver looks for a solution inside the specified interval and disregards any solutions it finds outside of the interval. When the solver finds the first solution inside the interval, it stops and does not look for other solutions. If you use MultiSolutions, the solver returns the solutions found outside of a specified interval.

    Note:   If you use the option MultiSolutions and do not specify any interval, the numeric solver returns only the first solution it finds.

With the MultiSolutions option, the solver also stops after it finds the first solution inside the specified interval. For example, find several numeric approximations for the following system:

eqs := [x*sin(10*x) = y^3, y^2 = exp(-2*x/3)]:
plot(x*sin(10*x) = y^3, y^2 = exp(-2*x/3))

Specify the interval where you want to search for the solutions. For example, consider the interval x = 0..1 that contains two solutions:

plot(x*sin(10*x) = y^3, y^2 = exp(-2*x/3), x = 0..1, y = 0..1)

Call the numeric::solve or numeric::fsolve command with the MultiSolutions option. Both solvers return one solution that belongs to the specified interval and one solution outside of the interval:

numeric::fsolve(eqs, [x = 0..1, y = 0..1], MultiSolutions);
numeric::solve(eqs, [x = 0..1, y = 0..1], MultiSolutions)

Specifying the interval that does not contain any solutions can help you find more approximations. In this case, the solver cannot find the solution inside the interval and continues searching. Before the solver quits, it can find many solutions outside the specified interval:

numeric::fsolve(eqs, [x = -10..0, y = 0..1], MultiSolutions)

Isolate Numeric Roots

If you know an interval containing the solutions of a polynomial equation, you can significantly speed up numeric approximations. To find one solution inside a specified interval, use the numeric::realroot command. The command returns real solutions and omits the complex ones:

numeric::realroot(1/4*x^4 + x^3 + x + 1 = 0, x = -5..0)

You also can specify an interval and search for all subintervals that can contain real roots. The numeric::realroots command returns a complete list of such subintervals:

numeric::realroots(1/4*x^4 + x^3 + x + 1 = 0, x = -5..0)

If the equation you solve is polynomial, each subinterval contains exactly one root. For nonpolynomial equations, numeric::realroots can return subintervals that do not contain any roots. numeric::realroots guarantees that the search interval does not contain any real roots outside the returned subintervals.

Solve Differential Equations and Systems

There is no general way to solve an arbitrary second- or higher-order ordinary differential equation (ODE). Often, such an equation does not have a symbolic solution. In some cases, the solution exists, but it can be presented only by special functions. For example, the solution of the following second-order ODE includes the error function erf:

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

Suppose, you do not need an exact symbolic solution, but you want to approximate the solution for several values of the parameter t. For numeric approximations of the solutions of ODEs, MuPAD provides two functions:

  • numeric::odesolve returns a numeric approximation of the solution at a particular point.

  • numeric::odesolve2 returns a function representing a numeric approximation of the solution.

Both functions accept either a first-order ODE or a system of first-order ODEs. To solve a higher-order equation, convert it to a system of the first-order equations. For example, represent the second-order ODE you solved symbolically as a system of two first-order equations: . The solution vector for this system is Y = [y, z].

Approximate at Particular Points

The function numeric::odesolve requires the system of first-order ODEs (dynamic system) to be represented by a procedure:

f := proc(t, Y) begin [Y[2], t*Y[2]] end_proc

The second parameter of numeric::odesolve is the range over which you want to solve an ODE. The third parameter is a list of initial conditions ( ). Approximate the solutions y(t) for t = 1, t = 3, and :

numeric::odesolve(f, 0..1, [0, 1/3]);
numeric::odesolve(f, 0..3, [0, 1/3]);
numeric::odesolve(f, 0..1/127, [0, 1/3])

MuPAD also offers an alternate way to generate parameters for the ODE numeric solvers:

  1. Define your initial value problem as a list or a set:

    IVP := {y''(t) = t*y'(t), y(0) = 0, y'(0) = 1/3}:
  2. Define a set of fields over which you want to get a solution:

    fields := [y(t), y'(t)]:
  3. Convert the initial value problem and the fields into a procedure acceptable by numeric::odesolve. The function numeric::ode2vectorfield generates the required procedure:

    [ODE, t0, Y0] := [numeric::ode2vectorfield(IVP, fields)]

Now call numeric::odesolve to approximate the solution at particular values of t:

numeric::odesolve(ODE, t0..1, Y0)

Represent Numeric Approximations as Functions

The function numeric::odesolve2 also requires the system of first-order ODEs (dynamic system) to be represented by a procedure. To generate parameters for numeric::odesolve2, use the following steps:

  1. Define your initial value problem as a list or a set:

    IVP := {y''(t) = t*y'(t), y(0) = 0, y'(0) = 1/3}:
  2. Define a set of fields over which you want to get a solution:

    fields := [y(t), y'(t)]:
  3. Convert the initial value problem and the fields into a procedure acceptable by numeric::odesolve2. The function numeric::ode2vectorfield generates the required procedure:

    ODE := numeric::ode2vectorfield(IVP, fields)

Now call numeric::odesolve2 to approximate the solution:

numApprox := numeric::odesolve2(ODE)

Using the function generated by numeric::odesolve2, find the numeric approximation at any point. For example, find the numeric solutions for the values t = 1, t = 3, and . You get the same results as with numeric::odesolve, but the syntax is shorter:

numApprox(1);
numApprox(3);
numApprox(1/127)

Plot the numeric solution using the function generated by numeric::odesolve2. The function numApprox returns a list [y(t), y'(t)]. When plotting the solution y(t), use brackets to extract the first entry of the solution list:

plotfunc2d(numApprox(t)[1], t = 0..3)

Use numeric::odesolve2 to find numeric approximations for the following system of ODEs:

IVP := {x'(t) = -y(t) + x(t)^2,
      y'(t) = 10*x(t) - y(t)^2,
            x(0) = 1, y(0) = 1}:
fields := [x(t), y(t)]:
ODE := numeric::ode2vectorfield(IVP, fields):
numApprox := numeric::odesolve2(ODE)

Plot the numeric solutions for x(t) and y(t) in one graph:

plotfunc2d(numApprox(t)[1], numApprox(t)[2], t = 0..20)

Use the plot::Curve2d plotting function to generate a parametric plot of the numeric solution:

curve := plot::Curve2d([numApprox(t)[1], numApprox(t)[2]],
                                               t = 0..20):
plot(curve)

Alternatively, use plot::Ode2d or plot::Ode3d.

Was this topic helpful?