Jonathan wrote:
> I'm trying to evaluate a nonlinear function to an accuraccy of 10^8. I
> have a function S1=1/(k^3x)^1/2 and I want to determine how many terms
> are needed to achieve the desired accuracy. To solve this requires an
> integration and an expansion using taylor series. My code is as follows:
>
> clc
> syms k x f=k^(3/2)/sqrt(1+x/k^3);
> y=taylor(f);
> z=int(y)
> a=subs(z,x,1) % subtitute x=1 into z
> a=1;
> while a<10^8;
> k=1;
> k=k+1;
> a=z+a;
> disp(k)
> disp(a)
> end
>
> How can I determine how many terms are required to get the desired
> accuracy? Thanks.
I am not going to answer that question directly, but I am going to point out
some items in what you have written.
Any given taylor series expansion has a fixed number of terms, and an error
order term. The evaluation of it at a particular point either is or is not
provably within the tolerance of the "real" solution of the untailored
formula. If it is not _provably_ within the tolerance of the real solution,
then there are two possibilities: A) that using a sufficiently higher order
taylor expansion will be provably within the desired tolerance; or B) that
even though the highest order term of the expansion may itself be smaller than
the desired tolerance, that the taylor series is in fact divergent and that
you cannot use a taylor expansion to calculate the function value to within
the desired tolerance.
It is not uncommon for people to write taylor expansion functions that iterate
until they find a term which is less than the tolerance value and then to stop
iterating  ignoring the possibility of divergence, or ignoring the issue
that a further term might be of a different sign and thus "widens" the current
window instead of narrowing it.
When you use a taylor series on a function, you do not need to integrate the
taylor series  not unless the function itself is defined by an integral, and
you are taking the taylor expansion of the integrand.
Now if you examine your code, you have y=taylor(f) . That computes a taylor
expansion to a fixed default order (5, according the documentation), around
the point where one of the variables is 0. You have not specified which
variable to operate on, so taylor() will call symvar, and will find that both
x and k are free variables in your formula  and will pick one of them as the
one to expand around. Which one will it pick? Well, symvar returns the
variables in alphabetical order, so probably taylor(f) is going to expand
around k=0 and treat x as a constant. Is that really what you want?
You never call taylor again, and you do not "manually" calculate taylor series
coefficients, so in view of what I have written above, you should now
understand that unless by chance the taylor expansion to the default order is
within the needed tolerance, then your routine in the form it is now will
never compute the answer correctly.
You have a line
a=subs(z,x,1) % subtitute x=1 into z
which has no obvious purpose. Why evaluate at 1 and not at 0 or 5 or 319?
Part of the problem here is that your problem statement is underspecified: you
have to achieve an accuracy of 1e8, but you have not specified the range over
which you must have that accuracy. The taylor expansion will be least accurate
at the extrema of the range that is furthest from the point of evaluation, so
if you are trying to prove something about the accuracy over the entire target
range, you should be picking the extrema. If your target range happens to be
[0,1] then you would already be doing that... except for the problem that
taylor() would likely have picked k=0 to expand around, not x=0 .
In the next statement of your code, you have the additional problem that you
completely overwrite the result of the substitution, with
a = 1. Note that a subs() operation does not affect the original expression:
it creates a new expression, which was returned in the variable 'a' in the
previous line.
On the next line, you start
while a<10^8
but you just assigned 1 to a, and 1 is not less than 1e8, so the while loop
will not execute at all. Probably you meant > instead of < . But even if so,
remember what I wrote above common pitfalls... including alternating signs.
And indeed, your S function *does* alternate signs in the taylor expansion
terms, so you should be considering at the very least using an abs() in the
test, and you should be thinking about how to deal with divergence and
cancellation. For example, there are taylor series for which each individual
term is greater than a particular value, but for which it can be shown that
adjacent terms of alternating sign contribute an increasingly negligible total
value.
Inside the while loop, in adjacent statements, you initialize k to 1 and then
you increment k by 1. What's the point of that?
As to your question, "How can I determine how many terms are required to get
the desired accuracy?": the answer can be quite difficult to find. A couple of
weeks ago, I was doing a taylor series of a function and showed that at even
500 decimal places and taylor order 999, the taylor expansion was convergent
and not very large, and yet the symbolic integrator was claiming that the real
answer was infinity. In order to understand that, I had to look at the limit
of the function as x approached infinity, and found that the function
decreased exponentially to a minimum of 10000 at infinity. But of course, a
minimum of 10000 an infinite number of times gives you an infinite sum (or
integral), so the _appearance_ of convergence was an illusion. So I tried
figuring out how many terms would be needed in order to sum to a value twice
that of the _apparent_ convergence point... and discovered that I couldn't
compute the number with the symbolic package I was using, which turned out not
to be able to process numbers with more than
((((((10^10)^10)^10)^10)^10)^10)^10 digits . The tail of the taylor expansion
was monotonically decreasing, but the tail did not converge to zero at
infinite order.
In order to determine the number of terms you will need to use, you need to
examine the behavior of the function with increasing taylor order, and
_somehow_ find an upper bound on the contribution of the neglected
coefficients, and establish that that contribution is less than your required
accuracy.
