Solve exponential and trig function using roots?

How can i solve an equation, such as, F(x) = 500 + 200*(E^X)(SIN(X)+COS(X)), using only the roots command? Is it possible?

[Torsten]

What makes you think you could use a root finder for polynomial functions for a function that is not a polynomial ?

Answer 1

[John D'Errico]

Simple answer - No. Roots is designed to solve only one function, polynomial root finding. It does that well. You cannot make code designed to do one thing, and make it suddenly do your bidding, thus to solve other unrelated problems. Software simply does not function that way. And just because "rootfinding" is a general term that describes the problem of finding the root(s) of a function, the MATLAB function roots is designed to solve the specific problem of POLYNOMIAL rootfinding.

Sometimes, rarely, but sometimes, you can use mathematics. And you need to be lucky. There are special circumstances where you could essentially subvert roots to do your bidding. I'll set up a simple problem where I know the problem has only one root, AND that root is sufficiently close to zero, because this approach will only work on some problems with roots near zero. And it needs to be a function that, while nonlinear, is well-behaved in the area of zero, thus smooth and sufficiently differentiable. In fact, I'll pick a problem where I don't know the answer in advance.

fun = @(x) x.*exp(x) - 0.1;

Clearly, this problem is purely nonlinear, not polynomial. And it should have one root that is near zero, although it might have other roots in the complex plane.

fplot(fun,[-1,1])
grid on
yline(0);

So that should clearly satisfy my requirements. It appears to have a root near zero. The function is well-behaved near zero. It looks like the root should be a little less than 0.1. If I had to guess from the plot, maybe around x = 0.08. Now, I need to convert this into a problem of polynomial rootfinding. The trick here is to compute a truncated Taylor series for fun. As long as the problem has small roots, the Taylor series will be convergent, and the function will be well approximated by a polynomial.

syms X
P = taylor(X*exp(X) - 0.1,'order',12)
P =
X^11/3628800 + X^10/362880 + X^9/40320 + X^8/5040 + X^7/720 + X^6/120 + X^5/24 + X^4/6 + X^3/2 + X^2 + X - 1/10

So it appears the poly in P has coefficients, that for x less than 1, will all be small enough to result in a number that is effectively near eps(0.1) in double precision arithmetic. That is my goal here. I use eps(0.1) because that is the constant term in the Taylor series. And, for small x, x^11/3628800 will be very near zero. The simple conclusion is that this polynomial approximation is entirely sufficient in double precision for numbers near 0.08 or so.

We are almost there now.

roots(sym2poly(P))
ans =
           3.37474552801313 +      5.62566697408421i
           3.37474552801313 -      5.62566697408421i
         0.0674210754470835 +      4.95972699491001i
         0.0674210754470835 -      4.95972699491001i
          -1.82325146643646 +      3.75537506402692i
          -1.82325146643646 -      3.75537506402692i
          -3.70930435148572 +     0.907368484145826i
          -3.70930435148572 -     0.907368484145826i
          -2.95524904911849 +      2.45786901638727i
          -2.95524904911849 -      2.45786901638727i
         0.0912765271608622 +                     0i

fun(xr(end))
ans =
0.0000000000000000000084424336087639885294056854381626

So xr(1) is indeed a root of the function, at least in double precision arithmetic.

fzero(fun,1)
ans =
        0.0912765271608623

fzero agrees with this root. Of course, fzero was wildly better at the problem, because fzero did it with no mathematical gyrations needed on my part.

syms X
solve(X*exp(X) - 0.1)
ans =
lambertw(0, 1/10)

vpa(ans)
ans =
0.09127652716086226429989572142318

And of course, the symbolic toolbox was able to see right through my little problem to recognize that the Lambert W special function can solve that problem directly.

As well, the other "roots" found of that polynomial are not sufficiently close to zero in magnitude for the truncated series to have been adequate. So those other polynomial roots are spurious, created by the artifact of our having converted this into a polynomial rootfinding problem.

fun(xr)
ans =
                      0.0000000000000000000084424336087639885294056854381626
 - 0.17346336892629345558570527669881 + 0.057898566460439199385505832481867i
 - 0.17346336892629345558570527669881 - 0.057898566460439199385505832481867i
  - 0.061555263633523215171627170205188 + 0.1964015536446014250671935389054i
  - 0.061555263633523215171627170205188 - 0.1964015536446014250671935389054i
    0.49002694146614155215508046713113 + 0.32619727076572122734746005345231i
    0.49002694146614155215508046713113 - 0.32619727076572122734746005345231i
      5.0618410695450202098990467551793 - 1.2290204905417056031607008133457i
      5.0618410695450202098990467551793 + 1.2290204905417056031607008133457i
      178.39451541924459314828751128591 - 69.836694148864202944930691695663i
      178.39451541924459314828751128591 + 69.836694148864202944930691695663i

So while we have "succeeded" in using roots, I would hardly suggest this is a good solution that can be always employed. I imagine I could use mathematics to find other ways to convert a problem that is fundamentally a nonlinear one, into a polynomial rootfinding one. For example:

fun2 = @(x) exp(x.^2 - 2*x) - 137;
fplot(fun2,[-2,2])
grid on
yline(0);

So here the polynomial trick won't be as happy, because the root is just a bit too large for a truncated polynomial to be sufficient. I might get a decent approximation to a root, perhaps a good starting value for a rootfinder, but not terribly accurate in double precision. So this time, I'll transform the problem into one that does have a solution, by taking logs. Here, the natural log. Thus, if we have:

exp(X^2 - 2*X) == 137

then

X^2 - 2*X - log(137) == 0

I just took the natural log of both sides. You should see that this problem is even perfectly well represented as a polynomial one. Roots should be perfect here, and if there are multiple roots, they will be accurate.

roots([1 -2 -log(137)])
ans =
          3.43310109239796
         -1.43310109239796

In fact, we should recognize that roots has found both roots of the problem, and that had I used fplot to examine a wider interval, it would have identified a second solution. At the same time however, these roots would not have benn viable roots to find using the Taylor series method I first used. At least, not without getting tricky.

fplot(fun2,[-2,4])
grid 
yline(0);

Again, the set of problems where such a solution is admissable is small. And you need to know the mathematics to know when you can transform your problem into one where roots will be admissable.

Hmm. Any other methods? The Taylor series trick is a nice one to occasionally gain quite good starting values for a root finder. And the log transformation is nice, when it will result in a valid polynomial.

But there are other such special cases. I can think of at least a few more. For example, what are the roots of the function:

fun3 = @(theta) sin(theta).^3 - 2*sin(2*theta) + 4*sin(theta) - 1.5
fplot(fun3,[-pi,pi])
grid
yline(0);

So there appear to be two real solutions, one near 1, and the other near 3. The trick here is a change of variables. So we will define u=sin(theta). Also, remember that sin(2*theta)=2*sin(theta)*cos(theta). (The classic double angle formula for the sine function.)

So our problem now transforms to

u^3 - 2*2*u*sqrt(1-u^2) + 4*u - 1.5 == 0

We will want later to check for any spurious solutions so introduced in this process. First, convert it to a classic polynomial.

Pu = expand((u^3 + 4*u - 1.5)^2 - 16*u^2*(1-u^2))
Pu =
u^6 + 24*u^4 - 3*u^3 - 12*u + 9/4

Solve, using roots.

uroots = roots(sym2poly(Pu))
uroots =
        -0.0520493014209754 +      4.89998585087286i
        -0.0520493014209754 -      4.89998585087286i
         -0.421885671187952 +     0.690657847776231i
         -0.421885671187952 -     0.690657847776231i
          0.759519625304405 +                     0i
          0.188350319913449 +                     0i

Transform back to theta, and test to see if any spurious roots were introduced.

theta = asin(uroots)
theta =
        -0.0104074545472298 +      2.29268598758042i
        -0.0104074545472298 -      2.29268598758042i
          -0.34695354346433 +      0.68063633762051i
          -0.34695354346433 -      0.68063633762051i
          0.862574308620704 +                     0i
          0.189482130786505 +                     0i

fun3(theta)
ans =
           4.08150263138723 -      196.016426690153i
           4.08150263138723 +      196.016426690153i
       1.99840144432528e-15 +  4.44089209850063e-16i
       1.99840144432528e-15 -  4.44089209850063e-16i
      -1.33226762955019e-15 +                     0i
          -1.47983366793158 +                     0i

So solutions 3, 4, amd 5 seem to have been true roots of the original problem, and only the 5th solution was a real root.

theta = asin(uroots(5))
theta =
         0.862574308620704

So, can you use roots to solve other problems that are not polynomial in nature? Rarely. And only if you know enough about the mathematics involved that you know when roots can indeed be viable.