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:
where is the mistake

Subject: where is the mistake

From: Nena Kabranski

Date: 26 Feb, 2011 21:16:07

Message: 1 of 21

function y=intxexp(n)
% y=intxexp(n)
% n is positive integer
% y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].
intxex=1;
y=1
for i=2:n
    i*intxex
    intxex=exp(1)-i*intxex;
    y=intxex
end

>>intxexp(17)
0.1043 correct!

>>intxexp(18)
0.8417 wrong!

>>intxexp(19)
-13.2742 wrong!
Thank you!

Subject: where is the mistake

From: Roger Stafford

Date: 26 Feb, 2011 22:10:21

Message: 2 of 21

"Nena Kabranski" wrote in message <ikbqin$4gr$1@fred.mathworks.com>...
> function y=intxexp(n)
> % y=intxexp(n)
> % n is positive integer
> % y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].
> intxex=1;
> y=1
> for i=2:n
> i*intxex
> intxex=exp(1)-i*intxex;
> y=intxex
> end
>
> >>intxexp(17)
> 0.1043 correct!
>
> >>intxexp(18)
> 0.8417 wrong!
>
> >>intxexp(19)
> -13.2742 wrong!
> Thank you!
- - - - - - -
  You are running into round-off problems with matlab's 53-bit double precision accuracy. The general expression for your p(n) is:

 (-1)^n * (!n * exp(1) - n!)

where !n is the notation for "subfactorial". These subfactorials have the famous sequence of values: 0, 1, 2, 9, 44, 265, 1854, etc. It is known that

 !n = floor(n!/exp(1))

so that !n * exp(1) - n! remains a small quantity.

Whatever the rounding errors are in the initial values of your p(n), they become greatly magnified by the time you get to n = 18 and 19 where neither !n nor n! can be expressed as exact integers. This is reflected in the fact that the above !n * exp(1) - n! expression becomes the small difference between two very large quantities.

Roger Stafford

Subject: where is the mistake

From: Nena Kabranski

Date: 27 Feb, 2011 04:59:14

Message: 3 of 21

Thank you very much for your response. It was very thorough and clear. I hope you teach somewhere so that more people benefit from your ability to explain.

Subject: where is the mistake

From: Roger Stafford

Date: 27 Feb, 2011 11:30:23

Message: 4 of 21

"Nena Kabranski" wrote in message <ikcln2$k52$1@fred.mathworks.com>...
> Thank you very much for your response. It was very thorough and clear. I hope you teach somewhere so that more people benefit from your ability to explain.
- - - - - - - - -
  It should be pointed out that one way around your round-off error problem is to compute your iteration backwards. This way the effect of initial errors will decrease rather than increase as the iteration proceeds. Of course the difficulty with that is not knowing what initial value to use to start the iteration. The remedy for that is to start a number of entries beyond the point where sufficiently accurate results are to be obtained and begin with an initial zero. It can be shown that by abandoning the top 18, the remaining p values will be accurate out to something like 14 or 15 decimal places. (Actually 18 may be overly cautious.) You can use the following code to get p values from p(1) to p(n) for an arbitrary n:

 p = zeros(n+18,1); % The zero at the top end initiates the iteration
 for k = n+17:-1:1
  p(k) = (exp(1)-p(k+1))/(k+1); % Your iteration going backwards
 end
 p = p(1:n); % Discard the top 18 entries

  With regard to teaching, many years ago I did give lectures in some college-level math courses but I have been retired for a long time now, so this newsgroup acts as a substitute to keep the little octogenarian "gray cells" active.

Roger Stafford

Subject: where is the mistake

From: Husam Aldahiyat

Date: 27 Feb, 2011 14:40:36

Message: 5 of 21

If you have the symbolic math toolbox then you can get around this like so:

 function y=intxexp(n)

 intxex=sym('1');
 y=sym('1');
 for i=2:n
     i*intxex
     intxex=exp(sym('1'))-i*intxex;
     y=intxex
 end

Subject: where is the mistake

From: John D'Errico

Date: 27 Feb, 2011 15:42:10

Message: 6 of 21

"Husam Aldahiyat" wrote in message <ikdnp4$4i$1@fred.mathworks.com>...
> If you have the symbolic math toolbox then you can get around this like so:
>
> function y=intxexp(n)
>
> intxex=sym('1');
> y=sym('1');
> for i=2:n
> i*intxex
> intxex=exp(sym('1'))-i*intxex;
> y=intxex
> end

I could not resist trying out my new toolbox in MATLAB.

FPF is a class that allows you to work in floating point
arithmetic with a user specified precision. So here, with
50 digits of precision in all computations...

intxex=fpf(1,50);
y=1;
n = 20;
e1 = fpf('e',50);
for i=2:n
  intxex=e1-i*intxex;
  y=intxex
end


y =
    0.7182818284590452353602874713526624977572470936999
y =
    0.5634363430819095292794250572946750044855058126002
y =
    0.4645364561314071182425872421739624798152238432991
y =
    0.3955995478020096441473512604828500986811278772044
y =
    0.3446845416469873704761799084555619056704798304735
y =
    0.3054900369301336420270281121637291580638882803854
y =
    0.2743615330179760991440625740428292332461408506167
y =
    0.2490280312972603430637243049671993985419794381496
y =
    0.2280015154864418047230444216806685123374527122039
y =
    0.210265158108185383406798832865308862045267259457
y =
    0.1950999311608206344787014769689561532140399802159
y =
    0.1819827233683769871371682707562325059747273508932
y =
    0.1705237013017674154399316807654074141110641811951
y =
    0.1604263089325340037613122598715512860912843757734
y =
    0.1514608855385011751792913134078419202966970813255
y =
    0.1434467743045252573123351434193498527133967111664
y =
    0.1362398909775906037382548898043651489161062927047
y =
    0.1297238998848237643334445650697246683512275323106
y =
    0.1238038307625699486913961699581691307326964474879


As expected, it agrees with Roger's reversed series.

John

Subject: where is the mistake

From: Nasser M. Abbasi

Date: 27 Feb, 2011 16:13:26

Message: 7 of 21

On 2/27/2011 7:42 AM, John D'Errico wrote:

>
> I could not resist trying out my new toolbox in MATLAB.
>
> FPF is a class that allows you to work in floating point
> arithmetic with a user specified precision. So here, with
> 50 digits of precision in all computations...
>
> intxex=fpf(1,50);
> y=1;
> n = 20;
> e1 = fpf('e',50);
> for i=2:n
> intxex=e1-i*intxex;
> y=intxex
> end
>
>
> y =
> 0.7182818284590452353602874713526624977572470936999
> y =
> 0.5634363430819095292794250572946750044855058126002
> y =
> 0.4645364561314071182425872421739624798152238432991
> y =
> 0.3955995478020096441473512604828500986811278772044
> y =
> 0.3446845416469873704761799084555619056704798304735
> y =
> 0.3054900369301336420270281121637291580638882803854
> y =
> 0.2743615330179760991440625740428292332461408506167
> y =
> 0.2490280312972603430637243049671993985419794381496
> y =
> 0.2280015154864418047230444216806685123374527122039
> y =
> 0.210265158108185383406798832865308862045267259457
> y =
> 0.1950999311608206344787014769689561532140399802159
> y =
> 0.1819827233683769871371682707562325059747273508932
> y =
> 0.1705237013017674154399316807654074141110641811951
> y =
> 0.1604263089325340037613122598715512860912843757734
> y =
> 0.1514608855385011751792913134078419202966970813255
> y =
> 0.1434467743045252573123351434193498527133967111664
> y =
> 0.1362398909775906037382548898043651489161062927047
> y =
> 0.1297238998848237643334445650697246683512275323106
> y =
> 0.1238038307625699486913961699581691307326964474879
>
>
> As expected, it agrees with Roger's reversed series.
>
> John

Fyi, here is the result from Mathematica using 50 digits
precision, it matches your results. Note, while comparing,
Mathematica has 50 digits after the decimal point, the above
output has 49 digits after the decimal.


------------------------
y=1;
For[i=2,i<=20,i++,
{
    y = Exp[1]-i*y;
    Print[N[y,50]]
}
]
------------------------

0.71828182845904523536028747135266249775724709369996
0.56343634308190952927942505729467500448550581260008
0.46453645613140711824258724217396247981522384329964
0.39559954780200964414735126048285009868112787720178
0.34468454164698737047617990845556190567047983048929
0.30549003693013364202702811216372915806388828027495
0.27436153301797609914406257404282923324614085150038
0.24902803129726034306372430496719939854197943019658
0.22800151548644180472304442168066851233745279173416
0.21026515810818538340679883286530886204526638462423
0.19509993116082063447870147696895615321405047820923
0.18198272336837698713716827075623250597459087697995
0.17052370130176741543993168076540741411297481598071
0.16042630893253400376131225987155128606262485398932
0.15146088553850117517929131340784192075524942987092
0.14344677430452525731233514341934984491800678589438
0.13623989097759060373825488980436528923312494760116
0.12972389988482376433344456506972200232787308927791
0.12380383076256994869139616995822245119978530814180

--Nasser

Subject: where is the mistake

From: Roger Stafford

Date: 27 Feb, 2011 20:52:27

Message: 8 of 21

"Nasser M. Abbasi" <nma@12000.org> wrote in message <ikdt79$bdd$1@speranza.aioe.org>...
> On 2/27/2011 7:42 AM, John D'Errico wrote:
> > .......
> > I could not resist trying out my new toolbox in MATLAB.
> >
> > FPF is a class that allows you to work in floating point
> > arithmetic with a user specified precision. So here, with
> > 50 digits of precision in all computations...
> >.......
> ......
> Fyi, here is the result from Mathematica using 50 digits
> precision, it matches your results. Note, while comparing,
> Mathematica has 50 digits after the decimal point, the above
> output has 49 digits after the decimal.
> ......
- - - - - - - - - -
  I used the direct formula p(n) = (-1)^n * (!n * exp(1) - n!) with the symbolic toolbox set at a precision of 100 to compute p(17:20) accurately for comparison purposes and noticed a curious phenomenon. By n = 20 John's value has lost about 20 decimal places of accuracy from the original 50, which one would more or less expect from the intrinsic nature of the iteration. However Nassar's answers apparently remain accurate to the full 50 digits which is very surprising since he is using the same iteration. How did you manage that using only 50 digits precision, Nassar? Possibly you are getting more precision than you realize with Mathematica.

  Here are our respective results for n = 20:

John 0.1238038307625699486913961699581691307326964474879
Nassar 0.12380383076256994869139616995822245119978530814180
Roger .123803830762569948691396169958222451199785308141796375457917

  The significance of this as far as John's results are concerned is that if we were to continue the same iteration up to, say, n = 100, doing everything with only 50 decimal digit precision, the results would surely have deteriorated catastrophically. Nassar's results remain a puzzle to me but presumably for a sufficiently large n his results would also begin to show a rapid decline in accuracy using this algorithm.

Roger Stafford

Subject: where is the mistake

From: John D'Errico

Date: 27 Feb, 2011 21:02:05

Message: 9 of 21

"Roger Stafford" wrote in message <ikedib$582$1@fred.mathworks.com>...
> "Nasser M. Abbasi" <nma@12000.org> wrote in message <ikdt79$bdd$1@speranza.aioe.org>...
> > On 2/27/2011 7:42 AM, John D'Errico wrote:
> > > .......
> > > I could not resist trying out my new toolbox in MATLAB.
> > >
> > > FPF is a class that allows you to work in floating point
> > > arithmetic with a user specified precision. So here, with
> > > 50 digits of precision in all computations...
> > >.......
> > ......
> > Fyi, here is the result from Mathematica using 50 digits
> > precision, it matches your results. Note, while comparing,
> > Mathematica has 50 digits after the decimal point, the above
> > output has 49 digits after the decimal.
> > ......
> - - - - - - - - - -
> I used the direct formula p(n) = (-1)^n * (!n * exp(1) - n!) with the symbolic toolbox set at a precision of 100 to compute p(17:20) accurately for comparison purposes and noticed a curious phenomenon. By n = 20 John's value has lost about 20 decimal places of accuracy from the original 50, which one would more or less expect from the intrinsic nature of the iteration. However Nassar's answers apparently remain accurate to the full 50 digits which is very surprising since he is using the same iteration. How did you manage that using only 50 digits precision, Nassar? Possibly you are getting more precision than you realize with Mathematica.
>
> Here are our respective results for n = 20:
>
> John 0.1238038307625699486913961699581691307326964474879
> Nassar 0.12380383076256994869139616995822245119978530814180
> Roger .123803830762569948691396169958222451199785308141796375457917
>
> The significance of this as far as John's results are concerned is that if we were to continue the same iteration up to, say, n = 100, doing everything with only 50 decimal digit precision, the results would surely have deteriorated catastrophically. Nassar's results remain a puzzle to me but presumably for a sufficiently large n his results would also begin to show a rapid decline in accuracy using this algorithm.
>
> Roger Stafford

Yes, but then nothing stops me from using as many
extra digits as I need.

I'm not saying this is a good solution, as it is not. It is
rarely a good idea to substitute the brute force of high
precision for the finesse of good numerical analysis.

John

Subject: where is the mistake

From: Nasser M. Abbasi

Date: 27 Feb, 2011 22:40:54

Message: 10 of 21

On Feb 27, 12:52 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Nasser M. Abbasi" <n...@12000.org> wrote in message <ikdt79$bd...@speranza.aioe.org>...> On 2/27/2011 7:42 AM, John D'Errico wrote:
> > > .......
> > > I could not resist trying out my new toolbox in MATLAB.
>
> > > FPF is a class that allows you to work in floating point
> > > arithmetic with a user specified precision. So here, with
> > > 50 digits of precision in all computations...
> > >.......
> > ......
> > Fyi, here is the result from Mathematica using 50 digits
> > precision, it matches your results. Note, while comparing,
> > Mathematica has 50 digits after the decimal point, the above
> > output has 49 digits after the decimal.
> > ......
>


> - - - - - - - - - -
>   I used the direct formula p(n) = (-1)^n * (!n * exp(1) - n!) with the symbolic toolbox set at a precision of 100 to compute p(17:20) accurately for comparison purposes and noticed a curious phenomenon.  By n = 20 John's value has lost about 20 decimal places of accuracy from the original 50, which one would more or less expect from the intrinsic nature of the iteration.  However Nassar's answers apparently remain accurate to the full 50 digits which is very surprising since he is using the same iteration.  How did you manage that using only 50 digits precision, Nassar?  Possibly you are getting more precision than you realize with Mathematica.
>
>   Here are our respective results for n = 20:
>
> John   0.1238038307625699486913961699581691307326964474879
> Nassar 0.12380383076256994869139616995822245119978530814180
> Roger   .123803830762569948691396169958222451199785308141796375457917
>
>   The significance of this as far as John's results are concerned is that if we were to continue the same iteration up to, say, n = 100, doing everything with only 50 decimal digit precision, the results would surely have deteriorated catastrophically.  Nassar's results remain a puzzle to me but presumably for a sufficiently large n his results would also begin to show a rapid decline in accuracy using this algorithm.
>
> Roger Stafford

Hello Roger,

I am no expert on this, but Mathematica by default uses variable
precision for floating points, which
keeps tracks of precision at each step as described in this page. This
might be different from what you might have been thinking of.

http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html

This behaviour can be changed if needed.

But one can see the result of fixing the maximum precision to say 55
is that at one point, Mathematica now complained after 6 iterations:

----------------------------------
$MaxPrecision=55;
y=1;
For[i=2,i<=20,i++,
{y=Exp[1]-i*y;
Print[N[y,50]]}
]

0.71828182845904523536028747135266249775724709369996
0.56343634308190952927942505729467500448550581260008
0.46453645613140711824258724217396247981522384329964
0.39559954780200964414735126048285009868112787720178
0.34468454164698737047617990845556190567047983048929
0.30549003693013364202702811216372915806388828027495

$MaxPrecision::prec: In increasing internal precision while attempting
to evaluate
...
the limit $MaxPrecision .... was reached.
Increasing the value of $MaxPrecision may help resolve the
uncertainty. >>

0.2743615330179760991440625740428292332461408515004


$MaxPrecision::prec: .....

0.249028031297260343063724304967199398541979430197

etc...

When using Mathematica, and if I want to work with just
machinePrecision computation, I add this line to the top of my script:

$MinPrecision = $MachinePrecision; $MaxPrecision = $MachinePrecision;

--Nasser

Subject: where is the mistake

From: Nena Kabranski

Date: 28 Feb, 2011 05:32:04

Message: 11 of 21

Thank you very much all for your help. I am an enthusiastic novice and try to employ everything that I’ve learned from you (thank you very much!) and I really succeeded in some other problems. I was able to solve similar problems using Roger’s and Husam’s strategies, but here is one that none of the strategies help:
(It is about the Bessel functions, defined by J(n)=(1/pi) *( integral of cos (sin(x)-n*x)dx ) for x from 0 to pi.)

Now here is the relationship: J(n+1) = 2*n*J(n) – J (n-1).
I know the first 2 terms, they are given: J(0) = 0.7651976865 and J(1)= 0.4400505857. Find the first 20 terms (Including J(0), so I suppose the first 21).

If I go by the most logical way:

n=20;
bessef=zeros(1,n)';
bessef(1)=0.7651976865;
bessef(2)=0.4400505857;
for i=3:n
     bessef(i)=2*(i-2)*bessef(i-1)-bessef(i-2);
end
bessef

 the error is on the 12 position.

When I tried Roger’s strategy (or at least what I thought it was):

n=20;
bessef=zeros(1,n)';
bessef(n-1)=2^(-52);
for i=2:n-1
     bessef(n-i)=2*(n-i+1)*bessef(n-i+1)-bessef(n-i+2);
end
bessef

I get numbers that are not even closed to the truth!
Thank you!

Subject: where is the mistake

From: Husam Aldahiyat

Date: 28 Feb, 2011 05:43:04

Message: 12 of 21

Try my strategy, or better yet, do this:

>> doc bessel

Subject: where is the mistake

From: Matt J

Date: 28 Feb, 2011 08:16:22

Message: 13 of 21

"Nena Kabranski" wrote in message <ikbqin$4gr$1@fred.mathworks.com>...
> function y=intxexp(n)
> % y=intxexp(n)
> % n is positive integer
> % y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].
> intxex=1;
> y=1
> for i=2:n
> i*intxex
> intxex=exp(1)-i*intxex;
> y=intxex
> end
=============

Yet another alternative is to expand exp*(x) into a Taylor series T(x) and integrate
T(x)*x^n term by term. This leads to a sum of positive terms which is very numerically stable. Using 21 terms, for example, you would do


factorials=factorial(0:21);

intexp=@(n) sum( 1./( factorials(2:end) + n*factorials(1:end-1) ) );

for n=1:20
 p(n,1)=intexp(n);
 
end

Here are the results, which agree pretty closely with what the others have proposed:


p =

   1.000000000000000
   0.718281828459045
   0.563436343081910
   0.464536456131407
   0.395599547802010
   0.344684541646987
   0.305490036930134
   0.274361533017976
   0.249028031297260
   0.228001515486442
   0.210265158108185
   0.195099931160821
   0.181982723368377
   0.170523701301767
   0.160426308932534
   0.151460885538501
   0.143446774304525
   0.136239890977591
   0.129723899884824
   0.123803830762570

Subject: where is the mistake

From: Nena Kabranski

Date: 28 Feb, 2011 23:28:07

Message: 14 of 21

"Husam Aldahiyat" wrote in message <ikfcl8$6tq$1@fred.mathworks.com>...
> Try my strategy, or better yet, do this:
>
> >> doc bessel

Husam, I used your method too. I don’t see anything wrong! For i=3:20 I get:
 bessef =
 1.0e+011 *
  0.000000000007652
  0.000000000004401
 -0.000000000025254 + 0.000000000008801i
  0.000000000079013 - 0.000000000085712i
 -0.000000000119376 + 0.000000000492074i
 -0.000000000585659 - 0.000000002121336i
  0.000000006704682 + 0.000000006821953i
 -0.000000039876975 - 0.000000011757111i
  0.000000176317442 - 0.000000039547458i
 -0.000000586297878 + 0.000000522581827i
  0.000001123710417 - 0.000003223375605i
  0.000002538207420 + 0.000014618341427i
 -0.000040513222953 - 0.000050173575264i
  0.000259861834920 + 0.000105049513722i
 -0.001209033144171 + 0.000149699190215i
  0.004276872361335 - 0.003121912562923i
 -0.009654631175322 + 0.020891695784145i
 -0.007441739228334 - 0.099754132924301i
  0.238929853937261 + 0.363241357456389i
 -1.674760391433489 - 0.875351589026733i

I don't want to work with the bessel function. I just want to find the mistake in this sequence!

Subject: where is the mistake

From: Roger Stafford

Date: 1 Mar, 2011 00:31:24

Message: 15 of 21

"Nena Kabranski" wrote in message <ikfc0k$p7b$1@fred.mathworks.com>...
> Thank you very much all for your help. I am an enthusiastic novice and try to employ everything that I’ve learned from you (thank you very much!) and I really succeeded in some other problems. I was able to solve similar problems using Roger’s and Husam’s strategies, but here is one that none of the strategies help:
> (It is about the Bessel functions, defined by J(n)=(1/pi) *( integral of cos (sin(x)-n*x)dx ) for x from 0 to pi.)
>
> Now here is the relationship: J(n+1) = 2*n*J(n) – J (n-1).
> I know the first 2 terms, they are given: J(0) = 0.7651976865 and J(1)= 0.4400505857. Find the first 20 terms (Including J(0), so I suppose the first 21).
>
> If I go by the most logical way:
>
> n=20;
> bessef=zeros(1,n)';
> bessef(1)=0.7651976865;
> bessef(2)=0.4400505857;
> for i=3:n
> bessef(i)=2*(i-2)*bessef(i-1)-bessef(i-2);
> end
> bessef
>
> the error is on the 12 position.
>
> When I tried Roger’s strategy (or at least what I thought it was):
>
> n=20;
> bessef=zeros(1,n)';
> bessef(n-1)=2^(-52);
> for i=2:n-1
> bessef(n-i)=2*(n-i+1)*bessef(n-i+1)-bessef(n-i+2);
> end
> bessef
>
> I get numbers that are not even closed to the truth!
> Thank you!
- - - - - - - - - - -
  Nena, in suggesting the "backwards" technique on the integral p-series you inquired about in the beginning of this thread, I did not wish to imply that this was a general method applicable to all or even many types of recurrence relations. That particular recurrence has the very special property that only a single value, namely one, can be used for the initial p(1) that will lead to a series that remains within bounds, and in fact converges to zero. With any other value no matter how small the difference from one, the series if computed accurately will always diverge unboundedly (and therefore certainly not be a valid answer to your integration problem.) That means it is an inherently unstable kind of computation if begun at n = 1 proceeding forward, given that computers are sure to make rounding errors. It is also the reason for being able to start with a high value of n, assume a
approximate zero value for p(n), and work back to significantly earlier values (at least 18 earlier I stated) and attain a very close approximation. That makes this an unusual recurrence relation indeed.

  With your Bessel function recurrence no such property holds. It is of a very different type. You would have to use entirely different techniques to do the computation accurately. That is why the computation of Bessel functions is a non-trivial task, not to be attained in such an easy manner.

Roger Stafford

Subject: where is the mistake

From: Nena Kabranski

Date: 1 Mar, 2011 02:26:05

Message: 16 of 21

I guess I am not trying to solve those problems… I am trying to understand in which cases the MatLab is making mistakes. I cannot see a problem here, because the numbers are close to each other and they are not extremely small. Or are they?
Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.
Could you help me in this analysis?

Subject: where is the mistake

From: Roger Stafford

Date: 1 Mar, 2011 04:43:05

Message: 17 of 21

"Nena Kabranski" wrote in message <ikhlft$3h1$1@fred.mathworks.com>...
> I guess I am not trying to solve those problems… I am trying to understand in which cases the MatLab is making mistakes. I cannot see a problem here, because the numbers are close to each other and they are not extremely small. Or are they?
> Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.
> Could you help me in this analysis?
- - - - - - - - - - -
  Nena, I am trying to get you to see that it isn't simply Matlab that is making the errors you have found. With the algorithm you have used, any computing system possessing the same 53-bit accuracy limitation would make the same kinds of errors. This is inherent in the very nature of digital computing. It is the algorithm being used that must be considered to be inappropriate.

  In the case of your working your way upward from bessef(1) and bessef(2), since you specified them to ten digits, let us suppose that bessef(2) is too high by .0000000001 (tenth digit one too high.) Then in

 bessef(3)=2*1*bessef(2)-bessef(1)

bessef(3) will accordingly be too high by .0000000002. Then pursuing this further we get

 bessef( 4) too high by .0000000007
 bessef( 5) too high by .0000000040
 bessef( 6) too high by .0000000313
 bessef( 7) too high by .0000003090
 ......
 bessef(12) too high by .2915473799

In this rapid build-up the 12th point would be off by almost .3 simply due to this one original error in the tenth digit of bessef(2). The results have become almost total garbage after only nine steps.

  This is a sign of a bad algorithm, not the fault of the computer doing the processing. The recurrence equality itself is mathematically valid but, taking into consideration the limitations of computing techniques, not the correct way to find distant Bessel values. In numerical analysis one must learn to use great care in selecting an appropriate algorithm in making the best use of existing computing capabilities. In this case the fault lies with trying to use this particular recurrence relation in this manner.

  Your second attempt to work your way downwards from bessef(20) and bessef(19) using the totally fictitious values of zero and 2^(-52) will accomplish nothing whatever, since you are dealing with a totally different kind of recurrence relation. Just because it worked for your earlier p values does not mean that it would work for this Bessel recurrence, and in fact it most certainly does not.

Roger Stafford

Subject: where is the mistake

From: Matt J

Date: 1 Mar, 2011 14:02:05

Message: 18 of 21

"Roger Stafford" wrote in message <ikhtgp$7e0$1@fred.mathworks.com>...
>
> Nena, I am trying to get you to see that it isn't simply Matlab that is making the errors you have found. ... It is the algorithm being used that must be considered to be inappropriate.
===============

@Nena,

...which, incidentally, is why MATLAB provides dedicated functions BESSELI, BESSELJ, etc... for people who want to compute Bessel functions: In implementing these functions, Matlab developers have undoubtedly taken the trouble to use smarter algorithms than these recursion methods that you are considering.

Subject: where is the mistake

From: Steven_Lord

Date: 1 Mar, 2011 15:35:28

Message: 19 of 21



"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message news:ikhtgp$7e0$1@fred.mathworks.com...
> "Nena Kabranski" wrote in message <ikhlft$3h1$1@fred.mathworks.com>...
>> I guess I am not trying to solve those problems… I am trying to
>> understand in which cases the MatLab is making mistakes. I cannot see a
>> problem here, because the numbers are close to each other and they are
>> not extremely small. Or are they?
>> Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.
>> Could you help me in this analysis?
> - - - - - - - - - - -
> Nena, I am trying to get you to see that it isn't simply Matlab that is
> making the errors you have found. With the algorithm you have used, any
> computing system possessing the same 53-bit accuracy limitation would make
> the same kinds of errors. This is inherent in the very nature of digital
> computing. It is the algorithm being used that must be considered to be
> inappropriate.

*snip Roger's excellent example*

In fact, there's an entire branch of mathematics dedicated to this type of
problem.

http://en.wikipedia.org/wiki/Numerical_analysis

You might also find this Cleve's Corner article interesting and informative:

http://www.mathworks.com/company/newsletters/news_notes/clevescorner/winter02_cleve.html

as well as this often-cited paper:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.102.244

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: where is the mistake

From: Nena Kabranski

Date: 2 Mar, 2011 03:01:23

Message: 20 of 21

I understand if you’re annoyed with my questions. I’ve been in this field for 2 weeks now. I promise when I become a professional I’ll help you the same way as you do. :)
Sincerely,
Nena

Roger, did I tell you that you’re awesome! You explain magnificently! You probably want to consider an instructor position for some of the colleges in your area. You’ll be an asset to any Math department.

Subject: where is the mistake

From: Mark Shore

Date: 23 Mar, 2011 22:53:02

Message: 21 of 21

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ikdrci$jch$1@fred.mathworks.com>...

>
> I could not resist trying out my new toolbox in MATLAB.
>
> FPF is a class that allows you to work in floating point
> arithmetic with a user specified precision. So here, with
> 50 digits of precision in all computations...
>
> intxex=fpf(1,50);
> y=1;
> n = 20;
> e1 = fpf('e',50);
> for i=2:n
> intxex=e1-i*intxex;
> y=intxex
> end
>

Just wondering what particular toolbox this is from, assuming it's been released... I didn't find the FPF function in VPI or John's other recent FEX submissions.

Tags for this Thread

No tags are associated with 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