Can the function "Taylor" process the input type "Function handle" directly?

Like what I said, every time I use the "taylor" function, I need to convert the function to a symbolic equation, which isn't very convenient because sometims these two things can't be converted, eg:

receive_power_theo = @(misalignment_integ) integral(@(theta) integ_func(theta,misalignment_integ),0,2.*pi)*N_UCA./2./pi;

Apparently, integ_func is a function handle. And this can't be converted because of the integral term.

Any answer helps.

Thanks in advance!

Answer 1

[John D'Errico]

Can taylor handle some completely general function handle? No. Why not?

A Taylor series is easy to compute, in theory. The coefficients of the terms are just derivatives, of increasing order for each power of x, so as long as you can symbolically differentiate what you give it, and do so for arbitrarily high orders of differentiation as needed by the number of terms in your Taylor series there is no problem.

https://en.wikipedia.org/wiki/Taylor_series

The problem arises when you want to give the function taylor some general function handle. In this case, you have written a function handle that internally uses a call to integral, which does an ADAPTIVE numerical integration. Taylor cannot solve that problem. Sorry. In general, Taylor NEEDS a function it can symbolically differentiate. Even then, not all functions have Taylor series.

[Torsten]

And what about a numerical Taylor expansion ? Are there functions available in MATLAB to accomplish this task, or does it have to programmed ?

[祥宇 崔]

Actually, I have the same question as Torsten does. Could give me some clues?

[John D'Errico]

You can compute derivatives, well estimate them. For example, my own derivest code can estimate the derivative of a function, or higher level derivatives. You can find derivest on the file exchange. It uses adaptive, robuust methods that not only estimate the desired derivative using a variety of approximations, but then it gets tricky using a variation of a higher order Richardson extrapolation to estimate the error in that prediction.

So, what the heck, lets try it! See what happens, and why and when there are issues.

fun = @exp;

Now we all know the Taylor series expansion. for the sin function. In theory, globally convergent, but things d o go to hell for large argument, taking a huge number of terms before you might get convergence, and by then, the odds are good that double precision arithmetic will fail you utterly.

We know the Taylor series for exp(x) is just...

E(x) = 1 + f'(0)*x + f''(0)*x^2/2 + f'''(0)*x^3/6 + f''''(0)*x^4/24 + f'''''(0)*x^5/120 + f''''''(0)*x^6/720 + ...

I'll form a Taylor polynomial, but since I don't know when things will go to hell, I'll flip it around from the standard order we see for a polynomial, with the constant term first. The zero'th order term is easy to generate.

P(1) = fun(0)
P =
     1

I'll also maintain a measure of the uncertainty in the coefficients as we go along. That zero'th order coefficient is known as well as we can compute a number using double preciion, so I'll just

err(1) = eps;

Next, for the first order term.

[P(2),err(2)] = derivest(fun,0);
P(2) = P(2)/factorial(1)
P =
                         1         0.999999999999994
err =
      2.22044604925031e-16      2.73498595451069e-14

Derivest had no problem there, though even so, we see there is just a little crumbling around the edges. That linear coefficient is not exactly 1. Pretty close though, but the finite difference approximations used are still not perfect, even though I did use some nifty mathematics to do the very best I could.

[P(3),err(3)] = derivest(fun,0,'derivativeorder',2);
P(3) = P(3)/factorial(1)
P(3) = P(3)/factorial(2)
P =
    1         0.999999999999994         0.499999999999947 

Still no problem, but they were the easy ones. What happens when we try to go further out?

[P(4),err(4)] = derivest(fun,0,'derivativeorder',3);
P(4) = P(4)/factorial(3)
P =
   1         0.999999999999994         0.499999999999947         0.166666666669408

[P(5),err(5)] = derivest(fun,0,'derivativeorder',4);
P(5) = P(5)/factorial(4)
P =
Columns 1 through 4
   1         0.999999999999994         0.499999999999947         0.166666666669408
Column 5
    0.0416666666549087
    
err
err =
Columns 1 through 4
    2.22044604925031e-16      2.73498595451069e-14      2.26360365289265e-12      2.45705518576479e-10
Column 5
    4.06542487336491e-09

The coefficients of the Taylor polynomial are slowly getting a little less easy to estimate. You can see those error estimates are slowly starting to grow in a relative sense. And derivest is surprisingly good at this. I thought it would be tripping way sooner.

[P(6),err(6)] = derivest(fun,0,'derivativeorder',5);
Error using derivest>check_params
DerivativeOrder must be scalar, one of [1 2 3 4].
Error in derivest (line 212)
par = check_params(par); 

But finally, derivest fails. I never implemented the code to go beyond that point. Could I have done so? I think so, since the finite difference schemes were automatically generated on the fly by derivest internally. Could you write the code? Again, you might be able to do so, but derivest is actually pretty tricky in how it works, and the trick to do a multi-level Richardson extrapolation was one I developed for derivest myself.

Anyway, the exponential is a well behaved function. The resulting Taylor series is not too terribly poor, and if we tried to evaluate that polynomial, compared to the true 5 term truncated Taylor series expansion, we would see little difference. The problem is, that truncated Tarlor sereis is only valid for a rather small domain.

Anyway, a Taylor series is something I would use only under extreme circumstances, for limited domains. And unless the first few terms give me what I need, I'd look for something better.

I'm not sure this helps any. While yes, you CAN generate a numerical Taylor series, I could probably have chosen a worse example with only a bit of effort, something that would have failed quickly.

[Torsten]

Thank you very much. I think this can help many MATLAB users.

Your code seems to be a kind of "Nseries" for MATLAB:

https://reference.wolfram.com/language/NumericalCalculus/ref/NSeries.html

[Paul]

@John D'Errico

Would the derivest result improve using VPA? Or is there some fundamental limitation in estimating the derivatives regardless of precision?

[John D'Errico]

Interesting question. I doubt it would work with VPA, because it was never written in that context. So some of the derivative approximations were actually written using linear algebra to generate them. In theory however, it could do so, and then the Richardson extrapolations could be quite high order, as long as the function is itself infintely differentiable. So, yes, it could work nicely then. It hs been a long time since I wrote that code, but I recall the general ideas pretty well.

Working in double precision, there is sort of a sweet spot that I needed to work with. For example, suppose you want to approximate a derivative. Derivest uses a sequence of function evals, approaching arbitrarily closely to the point in question. I'll use sin(x) as an example.

fun = @sin;

Now assume we want to compute the derivative of fun at x0==1. Generate a geometric sequence of points, increasing in distance away from x0.

format long g
x0 = 1;
xseq = x0 + logspace(-15,1,10)'

xseq = 10×1

1 1.00000000000006 1.00000000000359 1.00000000021544 1.0000000129155 1.00000077426368 1.00004641588834 1.00278255940221 1.16681005372001 11

The problem is, if fun is a black box, we don't really know much about the shape of fun. So it is not clear where we need to best sample the function to get a good estimate of the derivative.

fpapprox = (fun(xseq(2:end)) - fun(xseq(1)))./(xseq(2:end) - xseq(1))

fpapprox = 9×1

0.541509433962264 0.540327564894932 0.540302350077195 0.540302308791902 0.540301980167428 0.540282776862684 0.539130887896137 0.467779763398602 -0.18414611913586

Now, we know the derivative of fun at x==1, is just

cos(x0)

ans = 
          0.54030230586814

So in fact, the best approximation for the derivative of fun at x=x0 (for THIS particular low order approximation for the derivative) seems to arise when the increment dx is around 1e-8. The problem is that double precision is failing us. If a high precision arithmetic had been employed, then we would have no problems for an infinitessimally small dx.

In fact, had a higher precision arithmetic utility been used there, I could have avoided much of the crap I went through in derivest to find the best possible estimates.

[祥宇 崔]

Thank you so much for your patient and valuable answer!