## `numeric::fsolve` -- search for a numerical root of a system of equations

#### Introduction

`numeric::fsolve``(eqs, ..)` returns a numerical approximation of a solution of the system of equations `eqs`.

#### Call(s)

`numeric::fsolve(eq, x <, Options>)`
`numeric::fsolve(eq, x = a <, Options>)`
`numeric::fsolve(eq, x = a..b <, Options>)`
`numeric::fsolve(eqs, [x1, x2, ..] <, Options>)`
`numeric::fsolve(eqs, {x1, x2, ..} <, Options>)`
`numeric::fsolve(eqs, [x1 = a1, x2 = a2, ..] <, Options>)`
`numeric::fsolve(eqs, {x1 = a1, x2 = a2, ..} <, Options>)`
`numeric::fsolve(eqs, [x1 = a1..b1, x2 = a2..b2, ..] <, Options>)`
`numeric::fsolve(eqs, {x1 = a1..b1, x2 = a2..b2, ..} <, Options>)`

#### Parameters

`eq` | - | an arithmetical expression or an equation in one indeterminate `x`. An expression `eq` is interpreted as the equation `eq = 0`.
`x` | - | an identifier or an indexed identifier to be solved for.
`a` | - | real or complex numerical starting value for the internal search. Typically, a crude approximation of the solution.
`a..b` | - | a range of numerical values defining a search interval for the numerical root.
`eqs` | - | a list or a set of expressions or equations in several indeterminates `x1, x2, ..`. Expressions are interpreted as homogeneous equations.
`x1, x2, ..` | - | identifiers or indexed identifiers to be solved for.
`a1, a2, ..` | - | real or complex numerical starting values for the internal search. Typically, crude approximations of solution.
`a1..b1, a2..b2, ..` | - | ranges of numerical values defining search intervals for the numerical root.

#### Options

`RestrictedSearch` | - | makes `numeric::fsolve` return only numerical roots in the user-defined search range `x = a..b` and `[x1 = a1..b1, x2 = a2..b2, ..]`, respectively. This is the default search strategy, if a search range with real range parameters is specified for at least one of the unknowns.
`UnrestrictedSearch` | - | allows `numeric::fsolve` to find and return solutions outside the specified search range. With this option, the search range is only used to choose random starting points for the internal numerical search.
`MultiSolutions` | - | makes `numeric::fsolve` return all solutions found in the internal search
`Random` | - | With this option, several calls to `numeric::fsolve` with the same input parameters may produce different roots.

#### Returns

A single numerical root is returned as a list of equations `[x = value]` or `[x1 = value1, x2 = value2, ..]`, respectively. `FAIL` is returned, if no solution is found. With the option MultiSolutions, sequences of solutions may be returned.

#### Side Effects

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

#### Related Functions

`linsolve`, `numeric::linsolve`, `numeric::realroot`, `numeric::realroots`, `numeric::polyroots`, `numeric::polysysroots`, `numeric::solve`, `polylib::realroots`, `solve`

#### Details

- This is MuPAD's numerical solver for non-linear systems of equations. By default, it returns only one numerical solution.
- The equations must not contain symbolic objects other than the unknowns that cannot be converted to numerical values via `float`. Symbolic objects such as `PI` or `sqrt(2)` etc. are accepted. The same holds true for starting values and search ranges. Search ranges may contain +-`infinity`. Cf. example 2.
- `numeric::fsolve` implements a purely numerical Newton type root search with a working precision set by the environment variable `DIGITS`. Well separated simple roots should be exact within this precision. However, multiple roots or badly separated roots may be computed with a restricted precision. Cf. example 3.
- For systems of equations, the expressions defining the equations must have a symbolic derivative.
- Overdetermined systems (i.e., more equations than indeterminates) are not accepted. However, there may be more indeterminates than equations. Cf. example 4.
- Specifying indeterminates `[x1, x2, ..]` without starting values or search ranges is equivalent to the search ranges

               [x1 = -infinity..infinity, x2 = -infinity..infinity, ..].

  Note, however, that the user should assist `numeric::fsolve` by providing specific search ranges whenever possible!
- For real equations and real starting points or search ranges, the internal Newton iteration will usually produce real values, i.e., `numeric::fsolve` searches for real roots only (unless square roots, logarithms etc. happen to produce complex values from real input). Use complex starting points or search ranges to search for complex roots of real equations. Cf. example 5.
- Starting values and search ranges can be mixed. Cf. example 6.
- Search ranges should only be provided, if a solution is known to exist inside the search range. Otherwise, the search may take some time before `numeric::fsolve` gives up.
- Specification of a search range primarily means that starting points from this range are used for the internal Newton search. For sufficiently small search ranges enclosing a solution the search will usually pick out this solution. However, it may also happen that the Newton iteration drifts towards other solutions. With the default search strategy RestrictedSearch, only solutions from the search range are accepted, even if solutions outside the search range are found internally. With the search strategy UnrestrictedSearch, any solution outside the search range is accepted and returned. Cf. example 7.
- If starting values for all indeterminates are provided, then a single Newton iteration with these initial data is launched. It either leads to a solution or `numeric::fsolve` gives up and returns `FAIL`. The same holds true if search ranges `x = a..a` or `[x1 = a1..a1, x2 = a2..a2, ..]` of zero length are specified.

  The risc of failure is high when providing bad starting values! Starting values are appropriate only if a sufficiently good approximation of the solution is known! On the other hand, providing good starting values is the fastest way to a solution. Cf. example 8.

- If at least one of the indeterminates has a non-trivial search range, then `numeric::fsolve` uses several Newton iterations with different starting values from the search range. Cf. example 9. Search ranges in conjunction with the option UnrestrictedSearch provide a higher chance of detecting roots than (bad) starting values!

- User defined assumptions such as `assume(x > 0)` are not taken into account in the numerical search! Provide search ranges instead! Cf. example 2.

- Convergence may be slow for multiple roots. Furthermore, `numeric::fsolve` may fail to detect such roots.

- `setuserinfo(numeric::fsolve, 3)` provides detailed information on the internal search.
- Use `linsolve` or `numeric::linsolve` for systems of linear equations. Use `numeric::realroots`, if all real roots of a single non-polynomial real equation in a finite range are desired. Use `polylib::realroots`, if all real roots of a real univariate polynomial are desired. Use `numeric::polyroots`, if all real and complex roots of a univariate polynomial are desired. Use `numeric::solve`, if all roots of a multivariate polynomial system are desired.

#### Option: RestrictedSearch

- This is the default search strategy, so there is no need to specify this option explicitly. With this search strategy, only solutions matching search ranges `x = a..b` or `[x1 = a1..b1, x2 = a2..b2, ..]` with real range parameters are returned. A real starting point is regarded as the search range `[-infinity, infinity]`. Complex starting points or complex range parameters automatically switch to the search strategy UnrestrictedSearch.
- Once a root with components (r.1, r2, ..( is found, it is checked whether a.i<= r.i <= b.i is satisfied. If the root is not inside the search range, the search is continued. Note that solutions outside the search range may be found internally. These may be accessed with the option MultiSolutions. Cf. example 7.

#### Option: UnrestrictedSearch

- This option switches off the search strategy RestrictedSearch. With UnrestrictedSearch, `numeric::fsolve` stops its internal search whenever a root is found, even if the root is not inside the specified search range. Starting points for the internal Newton search are taken from the search range.

#### Option: MultiSolutions

- This option only has an effect when used with the default search strategy RestrictedSearch. A sequence of all roots found in the internal search is returned. Cf. example 7.

#### Option: Random

- With this option, random starting values are chosen for the internal search. Consequently, calling `numeric::fsolve` several times with the same parameters may lead to different solutions. This may be useful when several roots of one and the same equation or set of equations are desired.

#### Example 1

We compute roots of the sine function:

    >> numeric::fsolve(sin(x) = 0, x)

                                [x = -226.1946711]

With the option Random, several calls may result in different roots:

    >> numeric::fsolve(sin(x), x, Random)

                                 [x = 97.38937226]

    >> numeric::fsolve(sin(x), x, Random)

                                 [x = 53.40707511]

Particular solutions can be chosen by an appropriate starting point close to the wanted solution, or by a search interval:

    >> numeric::fsolve(sin(x), x = 3),
       numeric::fsolve(sin(x), x = -4..-3)

                       [x = 3.141592653], [x = -3.141592654]

The solutions found by `numeric::fsolve` can be used in `subs` and `assign` to substitute or assign the indeterminates:

    >> eqs := [x^2 = sin(y), y^2 = cos(x)]:

    >> solution := numeric::fsolve(eqs, [x, y])

                       [x = -0.8517004887, y = 0.8116062152]

    >> eval(subs(eqs, solution))

            [0.7253937224 = 0.7253937224, 0.6587046485 = 0.6587046485]

    >> assign(solution): x, y

                            -0.8517004887, 0.8116062152

    >> delete eqs, solution, x, y:

#### Example 2

We demonstrate the use of search ranges. The following system has solutions with positive and negative `x`. The solution with `x > 0` is obtained with the search interval `x = 0..infinity`:

    >> numeric::fsolve([x^2 = exp(x*y), x^2 = y^2],
                       [x = 0..infinity, y])

                       [x = 0.7530891649, y = -0.753089165]

    >> numeric::fsolve([x^2 = exp(x*y), x^2 = y^2],
                       [x = -infinity..0, y])

                       [x = -0.753089165, y = 0.7530891649]

#### Example 3

Multiple roots can only be computed with a restricted precision:

    >> numeric::fsolve(expand((x - 1/3)^5), x = 0.3)

                                [x = 0.3333929968]

#### Example 4

The following system of equations is degenerate and has a 1-parameter family of solutions. Each call to `numeric::fsolve` picks out one random solution:

    >> numeric::fsolve([x^2 - y^2, x^2 - y^2], [x, y], Random) $ i=1..3

          [x = -140.1698476, y = 140.1698476],

             [x = 34.70258251, y = 34.70258251],

             [x = -29.16650501, y = 29.16650501]

The equation may also be specified as an underdetermined system:

    >> numeric::fsolve([x^2 - y^2], [x, y])

                        [x = -140.1698476, y = 140.1698476]

#### Example 5

The following equation has no real solution. Consequently, the numerical search with real starting values fails:

    >> numeric::fsolve(sin(x) + cos(x)^2 = 3, x)

                                       FAIL

With a complex starting value, a solution is found:

    >> numeric::fsolve(sin(x) + cos(x)^2 = 3, x = I)

                        [x = 0.2972513613 + 1.128383965 I]

Also complex search ranges may be specified. In the following, the internal starting point is a random value on the line from `2 + I` to `3 + I`:

    >> numeric::fsolve(sin(x) + cos(x)^2 = 3, x = 2 + I..3 + I)

                         [x = 2.844341292 + 1.128383965 I]

#### Example 6

Starting values and search intervals can be mixed:

    >> numeric::fsolve([x^2 + y^2 = 1, y^2 + z^2 = 1, x^2 + z^2 = 1],
                       [x = 1, y = 0..10, z])

              [x = 0.7071067812, y = 0.7071067812, z = 0.7071067812]

#### Example 7

With UnrestrictedSearch, search intervals are only used for choosing starting values for the internal Newton search. The numerical iteration may drift towards a solution outside the search range:

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

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

                        [x = 1.232766202, y = -0.663038602]

With the default strategy RestrictedSearch, only solutions inside the search range are accepted:

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

                       [x = 0.9816416007, y = -0.7209295436]

In the last search, also the previous solution outside the search range was found. With the option MultiSolutions, `numeric::fsolve` returns a sequence of all solutions that were found in the internal search:

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

          [x = 0.9816416007, y = -0.7209295436],

             [x = 1.232766202, y = -0.663038602]

    >> delete eqs:

#### Example 8

Usually, most of the time is spent internally searching for some (crude) approximations of the root. If high precision roots are required, it is recommended to compute first approximations with moderate values of `DIGITS` and use them as starting values for a refined search:

    >> eq := exp(-x) = x:

    >> DIGITS := 10: firstApprox := numeric::fsolve(eq, x)

                                [x = 0.5671432904]

This output is suitable as input defining a starting value for x:

    >> DIGITS := 1000: numeric::fsolve(eq, firstApprox)

          [x = 0.5671432904097838729999686622103555497538...]

    >> delete eq, firstApprox, DIGITS:

#### Example 9

Specifying starting values for the indeterminates launches a single Newton iteration. This may fail, if the starting values are not sufficiently close to the solution:

    >> eq := [x*y = x + y - 4, x/y = x - y + 4]:

    >> numeric::fsolve(eq, [x = 1, y = 1])

                                       FAIL

If a search range is specified for at least one of the unknowns, then several Newton iterations with random starting values in the search range are used, until a solution is found or until `numeric::fsolve` gives up:

    >> numeric::fsolve(eq, [x = 1, y = 0..10])

                          [x = 4.026449604e-14, y = 4.0]

    >> delete eq:

#### Background

- Internally the set of equations f(x)=0 is solved by a modified Newton iteration x -> x - t*(f'(x))^(-1)*f(x) with some adaptively chosen stepsize t. For degenerate or ill-conditioned Jacobians f' a minimization strategy for <f,f> is implemented. For scalar real equations, `numeric::realroot` is used, if a real finite search range is specified.

#### Changes

- new function
- The functionality of the function `numeric::fsolve` of previous MuPAD versions was moved to `numeric::realroots`.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000