http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865
MATLAB Central Newsreader  Nonlinear functions need to find the accuracy
Feed for thread: Nonlinear functions need to find the accuracy
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Thu, 08 Apr 2010 19:17:25 +0000
Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734341
Jonathan
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:<br>
<br>
clc<br>
syms k x <br>
f=k^(3/2)/sqrt(1+x/k^3);<br>
y=taylor(f);<br>
z=int(y)<br>
a=subs(z,x,1) % subtitute x=1 into z<br>
a=1;<br>
while a<10^8;<br>
k=1;<br>
k=k+1;<br>
a=z+a;<br>
disp(k)<br>
disp(a)<br>
end<br>
<br>
How can I determine how many terms are required to get the desired accuracy? Thanks.

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

Thu, 08 Apr 2010 21:13:06 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734374
Jonathan
Thank you Walter for your response. As you can probably see, I am not very savy when it comes to programming. Perhaps a better explanation of the problem is needed. I know the desired answer, I just don't know how to get there in Matlab.<br>
<br>
I have two infinite series S1 and S2. S1=1/sqrt(k^3x) and S2=1/sqrt(k^3+x). Each are evaluated from k=1 to infinity. For the first part of the problem I want to evaluate S1 and S2 individually for any value of x. What I should find is that Matlab will not be able to compute these functions, and will hang until I manually stop the code. The second part of the problem, I need to substract S2 from S1 then do a summation from k=1 to infinity. Matlab can find a solution when this done, and it terns out to take about 150 terms to reach the desired accuracy. <br>
<br>
Hopefully this helps in defining what I am looking for. Sorry I don't have a great deal of experience in Matlab.

Fri, 09 Apr 2010 00:21:31 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734434
Walter Roberson
Jonathan wrote:<br>
<br>
> I have two infinite series S1 and S2. S1=1/sqrt(k^3x) and <br>
> S2=1/sqrt(k^3+x). Each are evaluated from k=1 to infinity. For the <br>
> first part of the problem I want to evaluate S1 and S2 individually for <br>
> any value of x. What I should find is that Matlab will not be able to <br>
> compute these functions, and will hang until I manually stop the code. <br>
> The second part of the problem, I need to substract S2 from S1 then do a <br>
> summation from k=1 to infinity. Matlab can find a solution when this <br>
> done, and it terns out to take about 150 terms to reach the desired <br>
> accuracy. <br>
<br>
I am still working on this question; the answers I am coming up with so far <br>
are that S1 needs well in excess of 100000 terms in order to have an accuracy <br>
better than 1e8 .<br>
<br>
In looking at your problem, I noticed an unstated assumption about the range <br>
of x. Suppose x = 1, then for your first term, k=1, you would have<br>
1/sqrt(1^31) which is 1/0 . This implies that x = 1 is not in the allowed <br>
range  but in your previous discussion, you were substituting x=1 into the <br>
taylor() of your function, which logically you should not have done.<br>
<br>
x has an infinite number of prohibited positive values, one for each k^3. <br>
Furthermore, if x > 1, then for all k < x^(1/3), you will have sqrt() of a <br>
negative number. I doubt you want that. I thus deduce that the valid range for <br>
x is < 1 .<br>
<br>
This logic leaves open the question of whether x can be negative. By examining <br>
S2 and by symmetry, we deduce that the answer is very likely no, and thus <br>
deduce that x is at most in [0, 1) .<br>
<br>
But can x be 0 itself? There is no logical reason why not in the expression of <br>
the sequences. However, when we go to expand the sequences and do integration, <br>
the x in the denominator sum will transform into something that includes a <br>
division by x, so if we attempt to evaluate that further processing at x=0, we <br>
will get division by 0. Really this is a failure of the symbolic integration <br>
routines to take into account the boundary conditions, but we have to work <br>
with what we have, not with what we wish we had. Thus, _effectively_ we have <br>
to exclude 0 from the valid range for x, and then specialcase it. The <br>
infinite sum for x=0 is Zeta(3/2)<br>
<br>
As x approaches 1, the infinite sum grows larger. If x = 1  1/10^N then the <br>
sum grows to _approximately_ sqrt(10^N) . As x can grow arbitrarily close to <br>
0, the infinite sum can grow arbitrarily large. At this point we have to ask <br>
whether your criteria is for accuracy of 1 part in 10^8, or if your criteria <br>
for accuracy is for 8 decimal places. If the latter, 8 decimal places, then <br>
the answer becomes that there is no fixed number of terms that can achieve the <br>
required accuracy over all of (0,1) : take x 1000 times closer to 1 and you <br>
require at least one more term to achieve the same number of decimal places of <br>
accuracy.<br>
<br>
<br>
So, at this point, before going forward, we need some clarification from you <br>
about the exact range of x, and a clarification as to what your accuracy <br>
requirement is.<br>
<br>
<br>
Oh yes, I should mention: considering the large number of terms and the small <br>
values that can be involved, my calculations show that double precision is not <br>
sufficient to be able to numerically achieve accuracy of 1E8 for S1  at <br>
least not in a straightforward manner.

Fri, 09 Apr 2010 00:46:22 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734436
Jonathan
X can be any value of my choosing, but for the sake of simplicity, I have chosen x=1. Mathematically, I know I can break the summation into two parts. The first part would be the function evaluated from k=1 to N1. The second part would be an integration of the function from N to infinity. To define accuracy better, I want the second part (the integration part) to be equal to 10^8. I've taken another crack at it, seperated into three parts. In my code, I manipulate the function. Here is my code thus far:<br>
<br>
clc<br>
syms z k<br>
%<br>
x=1;<br>
f=(1+z)^(1/2)<br>
g=taylor(f,3);<br>
h=subs(g,x/k^3,z);<br>
i=int(h*k^(3/2),k);<br>
j=abs(i);<br>
%<br>
k1=(1z)^(1/2)<br>
l=taylor(k1,3);<br>
m=subs(l,x/k^3,z);<br>
n=int(m*k^(3/2),k);<br>
o=abs(n);<br>
%<br>
p=kf;<br>
q=taylor(p,3);<br>
r=subs(q,z,x/k^3);<br>
s=int(r*k^(3/2),k);<br>
t=abs(s)<br>
%<br>
k=1;<br>
t=0;<br>
while s<(1/10^8); <br>
k=k+1;<br>
s=t+s;<br>
end<br>
disp(s)<br>
disp(k)

Fri, 09 Apr 2010 03:18:09 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734459
Roger Stafford
"Jonathan " <jauni6@hotmail.com> wrote in message <hpltcu$hfb$1@fred.mathworks.com>...<br>
> X can be any value of my choosing, but for the sake of simplicity, I have chosen x=1.<br>
<br>
I cannot resist commenting on one of the functions f(z) = (1z)^(1/2) that you are dealing with, (the k for the moment being taken out of the discussion for the sake of simplicity.) As you may have observed, its Taylor expansion about z = 0 in powers of z is:<br>
<br>
f(z) = 1 + 1/2*z + 3/8*z^2 + 5/16*z^3 + 35/128*z^4 + 63/256*z^5 + ...<br>
<br>
(The general rule for these coefficients can readily be derived from the fractional form of the binomial theorem.)<br>
<br>
Though your description was not very clear, I think it amounted to asking how far out this series one must go to achieve an accurate approximation to the original f(z). Also the same question was apparently asked of the integral of f(z) which would be<br>
<br>
2*(1z)^(1/2) + C<br>
<br>
and its power series expansion in powers of z.<br>
<br>
The answer to both questions very definitely depends on the value you select for z. You might expect some trouble when z is very nearly equal to 1 because f(z) approaches infinity there. But for what range of z values would this series converge? The answer might surprise you. It converges for any z less than 1 and greater than 1. The same holds true for its integral series even though no infinity occurs there at z = 1. The closer one gets to either of these limits the more terms one has to take to get a good approximation. Past the limits, both these series are divergent. No matter how far out you go, the sum gets farther and farther away from the function's true value.<br>
<br>
The explanation for this behavior is easiest to understand if you extend the definition of these functions to the complex plane. In the complex plane the range of convergence of power series like this is always a circle (though sometimes the circle might have an infinite radius) centered at z = 0 (assuming you have powers of z.) It can be shown that this circle always extends out to the first singularity on its perimeter in the function or in any of its derivatives. In our case that occurred at z = 1. Hence beyond the opposite side of the circle at z = 1 the series also must become divergent. In the complex plane for f(z) any z whose absolute value exceeds 1 also causes the series to diverge since these are also outside the circle. Again the same holds for the integral of f(z). Just because the derivative of this integral, namely f(z) itself, has a singularity at z = 1, then <br>
the integral's series also must diverge. This is true in spite of the fact that the functions are perfectly well defined beyond such a circle of convergence. <br>
<br>
Well, that's enough complex analysis. The lesson for you here is that it most definitely *does* depend on what values you assign to z as to how far out in the power series one needs to go for good approximation and indeed as to whether the series will converge at all.<br>
<br>
Roger Stafford

Fri, 09 Apr 2010 22:31:42 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734774
Walter Roberson
Jonathan wrote:<br>
> X can be any value of my choosing, but for the sake of simplicity, I <br>
> have chosen x=1.<br>
<br>
Then you are certain to get nonsense answers out. There is no polynomial, <br>
piecewise or continuous, that is going to be able to give you accuracy to <br>
within 1e=8 of a value that is undefined.<br>
<br>
1/sqrt(k^3x) when k = 1 and x = 1 is 1/sqrt(1^3  1) = 1/sqrt(1  1) = 1/0 . <br>
And since you are doing a summation, this result is going to pollute the <br>
entire rest of your sum.<br>
<br>
Unless you are willing to have your summation include complex terms, your S1 <br>
must restrict x to be strictly less than 1, and your S2 must restrict x to be <br>
strictly greater than 1. And when you switch over to the integral, you may <br>
have to exclude x = 0 as well and handle that case specially.<br>
<br>
Under the restriction that x is in (1, 1), int(S1,k=N..infinity) is<br>
<br>
1/3 * I * (4 * 3^(3/4) * (I / x^(1/3))^(1/2) * (1 / x^(1/3) / (3 + I * <br>
3^(1/2)))^(1/2) * (I / x^(1/3))^(1/2) * (N^3  x)^(1/2) * EllipticK(1/2 * <br>
3^(1/2)  1/2 * I) + 2 * I * 3^(3/4) * (I / x^(1/3))^(1/2) * (1 / x^(1/3) / <br>
(3 + I * 3^(1/2)))^(1/2) * (I / x^(1/3))^(1/2) * (N^3  x)^(1/2) *<br>
EllipticK((3 * I + 3^(1/2)) / (3 + I * 3^(1/2))) + 3 * EllipticF(1/6 * 2^(1/2) <br>
* ((3 * x^(1/3) + I * 3^(1/2) * x^(1/3) + 2 * I * 3^(1/2) * N) / <br>
x^(1/3))^(1/2) * 3^(1/2) , 1/2 * I + 1/2 * 3^(1/2)) * (I * (2 * N  I * <br>
3^(1/2) * x^(1/3) + x^(1/3)) / x^(1/3))^(1/2) * ((2 * x^(1/3) + 2 * N) / <br>
x^(1/3) / (3 + I * 3^(1/2)))^(1/2) * (I * (2 * N + x^(1/3) + I * 3^(1/2) * <br>
x^(1/3)) / x^(1/3))^(1/2) * (I / (3 + I * 3^(1/2)))^(1/2)) * x^(1/3) / (N^3  <br>
x)^(1/2) / (I / (3 + I * 3^(1/2)))^(1/2)<br>
<br>
Nice and simple, eh? And don't you love the way it tosses in the imaginary <br>
unit all over the place? So in order to produce a real answer, all of the <br>
imaginaries have to cancel out  and to make things even happier, the <br>
Elliptic* functions can produce complex results when presented with purely <br>
real input.<br>
<br>
Unfortunately, some testing with x=1/2 shows that the above integral is more <br>
or less bull excrement. For x=1/2, N=40000000000000001 suffices to reach the <br>
accuracy of 1/10^8, whereas the above integral with the same parameters gives <br>
a result of approximately 2.73 + 4.72I .<br>
<br>
Sigh. Time to talk to the manufacturer of my symbolic manipulation package <br>
again :(

Sat, 10 Apr 2010 00:02:36 +0000
Re: Nonlinear functions need to find the accuracy
http://www.mathworks.com/matlabcentral/newsreader/view_thread/278865#734787
Walter Roberson
Walter Roberson wrote:<br>
<br>
> Unfortunately, some testing with x=1/2 shows that the above integral is <br>
> more or less bull excrement. For x=1/2, N=40000000000000001 suffices to <br>
> reach the accuracy of 1/10^8, whereas the above integral with the same <br>
> parameters gives a result of approximately 2.73 + 4.72I .<br>
<br>
I found a workaround in my symbolic package, and using the workaround, I <br>
have found that within the range x > 1, x < 1, the solution for the integral <br>
of S1 to be 1/10^8 or less is N = 4 * 10^16 + 1 for positive x, and N = 4 * <br>
10^16 for negative x (a difference of exactly 1.) The answer is remarkably <br>
invariant for x in that range: changes in x over that range affect the 41'st <br>
decimal place of the answer. even though the number of working decimal places <br>
required goes up as x approaches 1 very closely.<br>
<br>
By symmetry on x, the limits are reversed for S2.<br>
<br>
For x = 0, 4 * 10^16 is the required number of terms (the integral is 2/sqrt(N) )