Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Nonlinear functions need to find the accuracy

Subject: Nonlinear functions need to find the accuracy

From: Jonathan

Date: 8 Apr, 2010 19:17:25

Message: 1 of 8

I'm trying to evaluate a nonlinear function to an accuraccy of 10^-8. I have a function S1=1/(k^3-x)^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.

Subject: Nonlinear functions need to find the accuracy

From: Walter Roberson

Date: 8 Apr, 2010 20:40:11

Message: 2 of 8

Jonathan wrote:
> I'm trying to evaluate a nonlinear function to an accuraccy of 10^-8. I
> have a function S1=1/(k^3-x)^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 un-tailored
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 1e-8, 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 1e-8, 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.

Subject: Nonlinear functions need to find the accuracy

From: Jonathan

Date: 8 Apr, 2010 21:13:06

Message: 3 of 8

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.

I have two infinite series S1 and S2. S1=1/sqrt(k^3-x) 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.

Hopefully this helps in defining what I am looking for. Sorry I don't have a great deal of experience in Matlab.

Subject: Nonlinear functions need to find the accuracy

From: Walter Roberson

Date: 9 Apr, 2010 00:21:31

Message: 4 of 8

Jonathan wrote:

> I have two infinite series S1 and S2. S1=1/sqrt(k^3-x) 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.

I am still working on this question; the answers I am coming up with so far
are that S1 needs well in excess of 100000 terms in order to have an accuracy
better than 1e-8 .

In looking at your problem, I noticed an unstated assumption about the range
of x. Suppose x = 1, then for your first term, k=1, you would have
1/sqrt(1^3-1) which is 1/0 . This implies that x = 1 is not in the allowed
range -- but in your previous discussion, you were substituting x=1 into the
taylor() of your function, which logically you should not have done.

x has an infinite number of prohibited positive values, one for each k^3.
Furthermore, if x > 1, then for all k < x^(1/3), you will have sqrt() of a
negative number. I doubt you want that. I thus deduce that the valid range for
x is < 1 .

This logic leaves open the question of whether x can be negative. By examining
S2 and by symmetry, we deduce that the answer is very likely no, and thus
deduce that x is at most in [0, 1) .

But can x be 0 itself? There is no logical reason why not in the expression of
the sequences. However, when we go to expand the sequences and do integration,
the x in the denominator sum will transform into something that includes a
division by x, so if we attempt to evaluate that further processing at x=0, we
will get division by 0. Really this is a failure of the symbolic integration
routines to take into account the boundary conditions, but we have to work
with what we have, not with what we wish we had. Thus, _effectively_ we have
to exclude 0 from the valid range for x, and then special-case it. The
infinite sum for x=0 is Zeta(3/2)

As x approaches 1, the infinite sum grows larger. If x = 1 - 1/10^N then the
sum grows to _approximately_ sqrt(10^N) . As x can grow arbitrarily close to
0, the infinite sum can grow arbitrarily large. At this point we have to ask
whether your criteria is for accuracy of 1 part in 10^8, or if your criteria
for accuracy is for 8 decimal places. If the latter, 8 decimal places, then
the answer becomes that there is no fixed number of terms that can achieve the
  required accuracy over all of (0,1) : take x 1000 times closer to 1 and you
require at least one more term to achieve the same number of decimal places of
accuracy.


So, at this point, before going forward, we need some clarification from you
about the exact range of x, and a clarification as to what your accuracy
requirement is.


Oh yes, I should mention: considering the large number of terms and the small
values that can be involved, my calculations show that double precision is not
sufficient to be able to numerically achieve accuracy of 1E-8 for S1 -- at
least not in a straight-forward manner.

Subject: Nonlinear functions need to find the accuracy

From: Jonathan

Date: 9 Apr, 2010 00:46:22

Message: 5 of 8

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 N-1. 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:

clc
syms z k
%
x=1;
f=(1+z)^(1/2)
g=taylor(f,3);
h=subs(g,x/k^3,z);
i=int(h*k^(-3/2),k);
j=abs(i);
%
k1=(1-z)^(1/2)
l=taylor(k1,3);
m=subs(l,x/k^3,z);
n=int(m*k^(-3/2),k);
o=abs(n);
%
p=k-f;
q=taylor(p,3);
r=subs(q,z,x/k^3);
s=int(r*k^(-3/2),k);
t=abs(s)
%
k=1;
t=0;
while s<(1/10^8);
    k=k+1;
    s=t+s;
end
disp(s)
disp(k)

Subject: Nonlinear functions need to find the accuracy

From: Roger Stafford

Date: 9 Apr, 2010 03:18:09

Message: 6 of 8

"Jonathan " <jauni6@hotmail.com> wrote in message <hpltcu$hfb$1@fred.mathworks.com>...
> X can be any value of my choosing, but for the sake of simplicity, I have chosen x=1.
-----------------
  I cannot resist commenting on one of the functions f(z) = (1-z)^(-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:

 f(z) = 1 + 1/2*z + 3/8*z^2 + 5/16*z^3 + 35/128*z^4 + 63/256*z^5 + ...

(The general rule for these coefficients can readily be derived from the fractional form of the binomial theorem.)

  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

 -2*(1-z)^(1/2) + C

 and its power series expansion in powers of z.

  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.

  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
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.

  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.

Roger Stafford

Subject: Nonlinear functions need to find the accuracy

From: Walter Roberson

Date: 9 Apr, 2010 22:31:42

Message: 7 of 8

Jonathan wrote:
> X can be any value of my choosing, but for the sake of simplicity, I
> have chosen x=1.

Then you are certain to get nonsense answers out. There is no polynomial,
piecewise or continuous, that is going to be able to give you accuracy to
within 1e=8 of a value that is undefined.

1/sqrt(k^3-x) when k = 1 and x = 1 is 1/sqrt(1^3 - 1) = 1/sqrt(1 - 1) = 1/0 .
And since you are doing a summation, this result is going to pollute the
entire rest of your sum.

Unless you are willing to have your summation include complex terms, your S1
must restrict x to be strictly less than 1, and your S2 must restrict x to be
strictly greater than -1. And when you switch over to the integral, you may
have to exclude x = 0 as well and handle that case specially.

Under the restriction that x is in (-1, 1), int(S1,k=N..infinity) is

-1/3 * I * (-4 * 3^(3/4) * (I / x^(1/3))^(1/2) * (-1 / x^(1/3) / (3 + I *
3^(1/2)))^(1/2) * (-I / x^(1/3))^(1/2) * (N^3 - x)^(1/2) * EllipticK(1/2 *
3^(1/2) - 1/2 * I) + 2 * I * 3^(3/4) * (I / x^(1/3))^(1/2) * (-1 / x^(1/3) /
(3 + I * 3^(1/2)))^(1/2) * (-I / x^(1/3))^(1/2) * (N^3 - x)^(1/2) *
EllipticK((3 * I + 3^(1/2)) / (3 + I * 3^(1/2))) + 3 * EllipticF(1/6 * 2^(1/2)
* (-(-3 * x^(1/3) + I * 3^(1/2) * x^(1/3) + 2 * I * 3^(1/2) * N) /
x^(1/3))^(1/2) * 3^(1/2) , 1/2 * I + 1/2 * 3^(1/2)) * (I * (2 * N - I *
3^(1/2) * x^(1/3) + x^(1/3)) / x^(1/3))^(1/2) * (-(-2 * x^(1/3) + 2 * N) /
x^(1/3) / (3 + I * 3^(1/2)))^(1/2) * (-I * (2 * N + x^(1/3) + I * 3^(1/2) *
x^(1/3)) / x^(1/3))^(1/2) * (I / (3 + I * 3^(1/2)))^(1/2)) * x^(1/3) / (N^3 -
x)^(1/2) / (I / (3 + I * 3^(1/2)))^(1/2)

Nice and simple, eh? And don't you love the way it tosses in the imaginary
unit all over the place? So in order to produce a real answer, all of the
imaginaries have to cancel out -- and to make things even happier, the
Elliptic* functions can produce complex results when presented with purely
real input.

Unfortunately, some testing with x=1/2 shows that the above integral is more
or less bull excrement. For x=1/2, N=40000000000000001 suffices to reach the
accuracy of 1/10^8, whereas the above integral with the same parameters gives
a result of approximately -2.73 + 4.72I .

Sigh. Time to talk to the manufacturer of my symbolic manipulation package
again :(

Subject: Nonlinear functions need to find the accuracy

From: Walter Roberson

Date: 10 Apr, 2010 00:02:36

Message: 8 of 8

Walter Roberson wrote:

> Unfortunately, some testing with x=1/2 shows that the above integral is
> more or less bull excrement. For x=1/2, N=40000000000000001 suffices to
> reach the accuracy of 1/10^8, whereas the above integral with the same
> parameters gives a result of approximately -2.73 + 4.72I .

I found a work-around in my symbolic package, and using the work-around, I
have found that within the range x > -1, x < 1, the solution for the integral
of S1 to be 1/10^8 or less is N = 4 * 10^16 + 1 for positive x, and N = 4 *
10^16 for negative x (a difference of exactly 1.) The answer is remarkably
invariant for x in that range: changes in x over that range affect the 41'st
decimal place of the answer. even though the number of working decimal places
required goes up as x approaches 1 very closely.

By symmetry on x, the limits are reversed for S2.

For x = 0, 4 * 10^16 is the required number of terms (the integral is 2/sqrt(N) )

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us