http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553
MATLAB Central Newsreader  where is the mistake
Feed for thread: where is the mistake
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

Sat, 26 Feb 2011 21:16:07 +0000
where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821852
Nena Kabranski
function y=intxexp(n)<br>
% y=intxexp(n)<br>
% n is positive integer<br>
% y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].<br>
intxex=1;<br>
y=1<br>
for i=2:n<br>
i*intxex<br>
intxex=exp(1)i*intxex;<br>
y=intxex<br>
end<br>
<br>
>>intxexp(17)<br>
0.1043 correct!<br>
<br>
>>intxexp(18)<br>
0.8417 wrong!<br>
<br>
>>intxexp(19)<br>
13.2742 wrong!<br>
Thank you!

Sat, 26 Feb 2011 22:10:21 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821865
Roger Stafford
"Nena Kabranski" wrote in message <ikbqin$4gr$1@fred.mathworks.com>...<br>
> function y=intxexp(n)<br>
> % y=intxexp(n)<br>
> % n is positive integer<br>
> % y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].<br>
> intxex=1;<br>
> y=1<br>
> for i=2:n<br>
> i*intxex<br>
> intxex=exp(1)i*intxex;<br>
> y=intxex<br>
> end<br>
> <br>
> >>intxexp(17)<br>
> 0.1043 correct!<br>
> <br>
> >>intxexp(18)<br>
> 0.8417 wrong!<br>
> <br>
> >>intxexp(19)<br>
> 13.2742 wrong!<br>
> Thank you!<br>
      <br>
You are running into roundoff problems with matlab's 53bit double precision accuracy. The general expression for your p(n) is:<br>
<br>
(1)^n * (!n * exp(1)  n!)<br>
<br>
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 <br>
<br>
!n = floor(n!/exp(1))<br>
<br>
so that !n * exp(1)  n! remains a small quantity.<br>
<br>
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.<br>
<br>
Roger Stafford

Sun, 27 Feb 2011 04:59:14 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821902
Nena Kabranski
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.

Sun, 27 Feb 2011 11:30:23 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821927
Roger Stafford
"Nena Kabranski" wrote in message <ikcln2$k52$1@fred.mathworks.com>...<br>
> 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.<br>
        <br>
It should be pointed out that one way around your roundoff 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:<br>
<br>
p = zeros(n+18,1); % The zero at the top end initiates the iteration<br>
for k = n+17:1:1<br>
p(k) = (exp(1)p(k+1))/(k+1); % Your iteration going backwards<br>
end<br>
p = p(1:n); % Discard the top 18 entries<br>
<br>
With regard to teaching, many years ago I did give lectures in some collegelevel 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.<br>
<br>
Roger Stafford

Sun, 27 Feb 2011 14:40:36 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821958
Husam Aldahiyat
If you have the symbolic math toolbox then you can get around this like so:<br>
<br>
function y=intxexp(n)<br>
<br>
intxex=sym('1');<br>
y=sym('1');<br>
for i=2:n<br>
i*intxex<br>
intxex=exp(sym('1'))i*intxex;<br>
y=intxex<br>
end

Sun, 27 Feb 2011 15:42:10 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821966
John D'Errico
"Husam Aldahiyat" wrote in message <ikdnp4$4i$1@fred.mathworks.com>...<br>
> If you have the symbolic math toolbox then you can get around this like so:<br>
> <br>
> function y=intxexp(n)<br>
> <br>
> intxex=sym('1');<br>
> y=sym('1');<br>
> for i=2:n<br>
> i*intxex<br>
> intxex=exp(sym('1'))i*intxex;<br>
> y=intxex<br>
> end<br>
<br>
I could not resist trying out my new toolbox in MATLAB.<br>
<br>
FPF is a class that allows you to work in floating point<br>
arithmetic with a user specified precision. So here, with<br>
50 digits of precision in all computations...<br>
<br>
intxex=fpf(1,50);<br>
y=1;<br>
n = 20;<br>
e1 = fpf('e',50);<br>
for i=2:n<br>
intxex=e1i*intxex;<br>
y=intxex<br>
end<br>
<br>
<br>
y =<br>
0.7182818284590452353602874713526624977572470936999<br>
y =<br>
0.5634363430819095292794250572946750044855058126002<br>
y =<br>
0.4645364561314071182425872421739624798152238432991<br>
y =<br>
0.3955995478020096441473512604828500986811278772044<br>
y =<br>
0.3446845416469873704761799084555619056704798304735<br>
y =<br>
0.3054900369301336420270281121637291580638882803854<br>
y =<br>
0.2743615330179760991440625740428292332461408506167<br>
y =<br>
0.2490280312972603430637243049671993985419794381496<br>
y =<br>
0.2280015154864418047230444216806685123374527122039<br>
y =<br>
0.210265158108185383406798832865308862045267259457<br>
y =<br>
0.1950999311608206344787014769689561532140399802159<br>
y =<br>
0.1819827233683769871371682707562325059747273508932<br>
y =<br>
0.1705237013017674154399316807654074141110641811951<br>
y =<br>
0.1604263089325340037613122598715512860912843757734<br>
y =<br>
0.1514608855385011751792913134078419202966970813255<br>
y =<br>
0.1434467743045252573123351434193498527133967111664<br>
y =<br>
0.1362398909775906037382548898043651489161062927047<br>
y =<br>
0.1297238998848237643334445650697246683512275323106<br>
y =<br>
0.1238038307625699486913961699581691307326964474879<br>
<br>
<br>
As expected, it agrees with Roger's reversed series.<br>
<br>
John

Sun, 27 Feb 2011 16:13:26 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#821969
Nasser M. Abbasi
On 2/27/2011 7:42 AM, John D'Errico wrote:<br>
<br>
><br>
> I could not resist trying out my new toolbox in MATLAB.<br>
><br>
> FPF is a class that allows you to work in floating point<br>
> arithmetic with a user specified precision. So here, with<br>
> 50 digits of precision in all computations...<br>
><br>
> intxex=fpf(1,50);<br>
> y=1;<br>
> n = 20;<br>
> e1 = fpf('e',50);<br>
> for i=2:n<br>
> intxex=e1i*intxex;<br>
> y=intxex<br>
> end<br>
><br>
><br>
> y =<br>
> 0.7182818284590452353602874713526624977572470936999<br>
> y =<br>
> 0.5634363430819095292794250572946750044855058126002<br>
> y =<br>
> 0.4645364561314071182425872421739624798152238432991<br>
> y =<br>
> 0.3955995478020096441473512604828500986811278772044<br>
> y =<br>
> 0.3446845416469873704761799084555619056704798304735<br>
> y =<br>
> 0.3054900369301336420270281121637291580638882803854<br>
> y =<br>
> 0.2743615330179760991440625740428292332461408506167<br>
> y =<br>
> 0.2490280312972603430637243049671993985419794381496<br>
> y =<br>
> 0.2280015154864418047230444216806685123374527122039<br>
> y =<br>
> 0.210265158108185383406798832865308862045267259457<br>
> y =<br>
> 0.1950999311608206344787014769689561532140399802159<br>
> y =<br>
> 0.1819827233683769871371682707562325059747273508932<br>
> y =<br>
> 0.1705237013017674154399316807654074141110641811951<br>
> y =<br>
> 0.1604263089325340037613122598715512860912843757734<br>
> y =<br>
> 0.1514608855385011751792913134078419202966970813255<br>
> y =<br>
> 0.1434467743045252573123351434193498527133967111664<br>
> y =<br>
> 0.1362398909775906037382548898043651489161062927047<br>
> y =<br>
> 0.1297238998848237643334445650697246683512275323106<br>
> y =<br>
> 0.1238038307625699486913961699581691307326964474879<br>
><br>
><br>
> As expected, it agrees with Roger's reversed series.<br>
><br>
> John<br>
<br>
Fyi, here is the result from Mathematica using 50 digits<br>
precision, it matches your results. Note, while comparing,<br>
Mathematica has 50 digits after the decimal point, the above<br>
output has 49 digits after the decimal.<br>
<br>
<br>
<br>
y=1;<br>
For[i=2,i<=20,i++,<br>
{<br>
y = Exp[1]i*y;<br>
Print[N[y,50]]<br>
}<br>
]<br>
<br>
<br>
0.71828182845904523536028747135266249775724709369996<br>
0.56343634308190952927942505729467500448550581260008<br>
0.46453645613140711824258724217396247981522384329964<br>
0.39559954780200964414735126048285009868112787720178<br>
0.34468454164698737047617990845556190567047983048929<br>
0.30549003693013364202702811216372915806388828027495<br>
0.27436153301797609914406257404282923324614085150038<br>
0.24902803129726034306372430496719939854197943019658<br>
0.22800151548644180472304442168066851233745279173416<br>
0.21026515810818538340679883286530886204526638462423<br>
0.19509993116082063447870147696895615321405047820923<br>
0.18198272336837698713716827075623250597459087697995<br>
0.17052370130176741543993168076540741411297481598071<br>
0.16042630893253400376131225987155128606262485398932<br>
0.15146088553850117517929131340784192075524942987092<br>
0.14344677430452525731233514341934984491800678589438<br>
0.13623989097759060373825488980436528923312494760116<br>
0.12972389988482376433344456506972200232787308927791<br>
0.12380383076256994869139616995822245119978530814180<br>
<br>
Nasser

Sun, 27 Feb 2011 20:52:27 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822000
Roger Stafford
"Nasser M. Abbasi" <nma@12000.org> wrote in message <ikdt79$bdd$1@speranza.aioe.org>...<br>
> On 2/27/2011 7:42 AM, John D'Errico wrote:<br>
> > .......<br>
> > I could not resist trying out my new toolbox in MATLAB.<br>
> ><br>
> > FPF is a class that allows you to work in floating point<br>
> > arithmetic with a user specified precision. So here, with<br>
> > 50 digits of precision in all computations...<br>
> >.......<br>
> ......<br>
> Fyi, here is the result from Mathematica using 50 digits<br>
> precision, it matches your results. Note, while comparing,<br>
> Mathematica has 50 digits after the decimal point, the above<br>
> output has 49 digits after the decimal.<br>
> ......<br>
         <br>
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.<br>
<br>
Here are our respective results for n = 20:<br>
<br>
John 0.1238038307625699486913961699581691307326964474879<br>
Nassar 0.12380383076256994869139616995822245119978530814180<br>
Roger .123803830762569948691396169958222451199785308141796375457917<br>
<br>
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.<br>
<br>
Roger Stafford

Sun, 27 Feb 2011 21:02:05 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822001
John D'Errico
"Roger Stafford" wrote in message <ikedib$582$1@fred.mathworks.com>...<br>
> "Nasser M. Abbasi" <nma@12000.org> wrote in message <ikdt79$bdd$1@speranza.aioe.org>...<br>
> > On 2/27/2011 7:42 AM, John D'Errico wrote:<br>
> > > .......<br>
> > > I could not resist trying out my new toolbox in MATLAB.<br>
> > ><br>
> > > FPF is a class that allows you to work in floating point<br>
> > > arithmetic with a user specified precision. So here, with<br>
> > > 50 digits of precision in all computations...<br>
> > >.......<br>
> > ......<br>
> > Fyi, here is the result from Mathematica using 50 digits<br>
> > precision, it matches your results. Note, while comparing,<br>
> > Mathematica has 50 digits after the decimal point, the above<br>
> > output has 49 digits after the decimal.<br>
> > ......<br>
>          <br>
> 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.<br>
> <br>
> Here are our respective results for n = 20:<br>
> <br>
> John 0.1238038307625699486913961699581691307326964474879<br>
> Nassar 0.12380383076256994869139616995822245119978530814180<br>
> Roger .123803830762569948691396169958222451199785308141796375457917<br>
> <br>
> 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.<br>
> <br>
> Roger Stafford<br>
<br>
Yes, but then nothing stops me from using as many<br>
extra digits as I need.<br>
<br>
I'm not saying this is a good solution, as it is not. It is<br>
rarely a good idea to substitute the brute force of high<br>
precision for the finesse of good numerical analysis.<br>
<br>
John

Sun, 27 Feb 2011 22:40:54 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822009
Nasser M. Abbasi
On Feb 27, 12:52 pm, "Roger Stafford"<br>
<ellieandrogerxy...@mindspring.com.invalid> wrote:<br>
> "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:<br>
> > > .......<br>
> > > I could not resist trying out my new toolbox in MATLAB.<br>
><br>
> > > FPF is a class that allows you to work in floating point<br>
> > > arithmetic with a user specified precision. So here, with<br>
> > > 50 digits of precision in all computations...<br>
> > >.......<br>
> > ......<br>
> > Fyi, here is the result from Mathematica using 50 digits<br>
> > precision, it matches your results. Note, while comparing,<br>
> > Mathematica has 50 digits after the decimal point, the above<br>
> > output has 49 digits after the decimal.<br>
> > ......<br>
><br>
<br>
<br>
>          <br>
> 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.<br>
><br>
> Here are our respective results for n = 20:<br>
><br>
> John 0.1238038307625699486913961699581691307326964474879<br>
> Nassar 0.12380383076256994869139616995822245119978530814180<br>
> Roger .123803830762569948691396169958222451199785308141796375457917<br>
><br>
> 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.<br>
><br>
> Roger Stafford<br>
<br>
Hello Roger,<br>
<br>
I am no expert on this, but Mathematica by default uses variable<br>
precision for floating points, which<br>
keeps tracks of precision at each step as described in this page. This<br>
might be different from what you might have been thinking of.<br>
<br>
<a href="http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html">http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html</a><br>
<br>
This behaviour can be changed if needed.<br>
<br>
But one can see the result of fixing the maximum precision to say 55<br>
is that at one point, Mathematica now complained after 6 iterations:<br>
<br>
<br>
$MaxPrecision=55;<br>
y=1;<br>
For[i=2,i<=20,i++,<br>
{y=Exp[1]i*y;<br>
Print[N[y,50]]}<br>
]<br>
<br>
0.71828182845904523536028747135266249775724709369996<br>
0.56343634308190952927942505729467500448550581260008<br>
0.46453645613140711824258724217396247981522384329964<br>
0.39559954780200964414735126048285009868112787720178<br>
0.34468454164698737047617990845556190567047983048929<br>
0.30549003693013364202702811216372915806388828027495<br>
<br>
$MaxPrecision::prec: In increasing internal precision while attempting<br>
to evaluate<br>
...<br>
the limit $MaxPrecision .... was reached.<br>
Increasing the value of $MaxPrecision may help resolve the<br>
uncertainty. >><br>
<br>
0.2743615330179760991440625740428292332461408515004<br>
<br>
<br>
$MaxPrecision::prec: .....<br>
<br>
0.249028031297260343063724304967199398541979430197<br>
<br>
etc...<br>
<br>
When using Mathematica, and if I want to work with just<br>
machinePrecision computation, I add this line to the top of my script:<br>
<br>
$MinPrecision = $MachinePrecision; $MaxPrecision = $MachinePrecision;<br>
<br>
Nasser

Mon, 28 Feb 2011 05:32:04 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822056
Nena Kabranski
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: <br>
(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.)<br>
<br>
Now here is the relationship: J(n+1) = 2*n*J(n) – J (n1). <br>
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).<br>
<br>
If I go by the most logical way:<br>
<br>
n=20;<br>
bessef=zeros(1,n)';<br>
bessef(1)=0.7651976865;<br>
bessef(2)=0.4400505857;<br>
for i=3:n<br>
bessef(i)=2*(i2)*bessef(i1)bessef(i2);<br>
end<br>
bessef<br>
<br>
the error is on the 12 position. <br>
<br>
When I tried Roger’s strategy (or at least what I thought it was):<br>
<br>
n=20;<br>
bessef=zeros(1,n)';<br>
bessef(n1)=2^(52);<br>
for i=2:n1<br>
bessef(ni)=2*(ni+1)*bessef(ni+1)bessef(ni+2);<br>
end<br>
bessef<br>
<br>
I get numbers that are not even closed to the truth!<br>
Thank you!

Mon, 28 Feb 2011 05:43:04 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822059
Husam Aldahiyat
Try my strategy, or better yet, do this:<br>
<br>
>> doc bessel

Mon, 28 Feb 2011 08:16:22 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822077
Matt J
"Nena Kabranski" wrote in message <ikbqin$4gr$1@fred.mathworks.com>...<br>
> function y=intxexp(n)<br>
> % y=intxexp(n)<br>
> % n is positive integer<br>
> % y calculates p(n) the integrals of (x^n)*exp(x) for x is [0,1].<br>
> intxex=1;<br>
> y=1<br>
> for i=2:n<br>
> i*intxex<br>
> intxex=exp(1)i*intxex;<br>
> y=intxex<br>
> end<br>
=============<br>
<br>
Yet another alternative is to expand exp*(x) into a Taylor series T(x) and integrate <br>
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<br>
<br>
<br>
factorials=factorial(0:21);<br>
<br>
intexp=@(n) sum( 1./( factorials(2:end) + n*factorials(1:end1) ) );<br>
<br>
for n=1:20<br>
p(n,1)=intexp(n);<br>
<br>
end<br>
<br>
Here are the results, which agree pretty closely with what the others have proposed:<br>
<br>
<br>
p =<br>
<br>
1.000000000000000<br>
0.718281828459045<br>
0.563436343081910<br>
0.464536456131407<br>
0.395599547802010<br>
0.344684541646987<br>
0.305490036930134<br>
0.274361533017976<br>
0.249028031297260<br>
0.228001515486442<br>
0.210265158108185<br>
0.195099931160821<br>
0.181982723368377<br>
0.170523701301767<br>
0.160426308932534<br>
0.151460885538501<br>
0.143446774304525<br>
0.136239890977591<br>
0.129723899884824<br>
0.123803830762570

Mon, 28 Feb 2011 23:28:07 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822301
Nena Kabranski
"Husam Aldahiyat" wrote in message <ikfcl8$6tq$1@fred.mathworks.com>...<br>
> Try my strategy, or better yet, do this:<br>
> <br>
> >> doc bessel<br>
<br>
Husam, I used your method too. I don’t see anything wrong! For i=3:20 I get:<br>
bessef =<br>
1.0e+011 *<br>
0.000000000007652 <br>
0.000000000004401 <br>
0.000000000025254 + 0.000000000008801i<br>
0.000000000079013  0.000000000085712i<br>
0.000000000119376 + 0.000000000492074i<br>
0.000000000585659  0.000000002121336i<br>
0.000000006704682 + 0.000000006821953i<br>
0.000000039876975  0.000000011757111i<br>
0.000000176317442  0.000000039547458i<br>
0.000000586297878 + 0.000000522581827i<br>
0.000001123710417  0.000003223375605i<br>
0.000002538207420 + 0.000014618341427i<br>
0.000040513222953  0.000050173575264i<br>
0.000259861834920 + 0.000105049513722i<br>
0.001209033144171 + 0.000149699190215i<br>
0.004276872361335  0.003121912562923i<br>
0.009654631175322 + 0.020891695784145i<br>
0.007441739228334  0.099754132924301i<br>
0.238929853937261 + 0.363241357456389i<br>
1.674760391433489  0.875351589026733i<br>
<br>
I don't want to work with the bessel function. I just want to find the mistake in this sequence!

Tue, 01 Mar 2011 00:31:24 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822310
Roger Stafford
"Nena Kabranski" wrote in message <ikfc0k$p7b$1@fred.mathworks.com>...<br>
> 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: <br>
> (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.)<br>
> <br>
> Now here is the relationship: J(n+1) = 2*n*J(n) – J (n1). <br>
> 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).<br>
> <br>
> If I go by the most logical way:<br>
> <br>
> n=20;<br>
> bessef=zeros(1,n)';<br>
> bessef(1)=0.7651976865;<br>
> bessef(2)=0.4400505857;<br>
> for i=3:n<br>
> bessef(i)=2*(i2)*bessef(i1)bessef(i2);<br>
> end<br>
> bessef<br>
> <br>
> the error is on the 12 position. <br>
> <br>
> When I tried Roger’s strategy (or at least what I thought it was):<br>
> <br>
> n=20;<br>
> bessef=zeros(1,n)';<br>
> bessef(n1)=2^(52);<br>
> for i=2:n1<br>
> bessef(ni)=2*(ni+1)*bessef(ni+1)bessef(ni+2);<br>
> end<br>
> bessef<br>
> <br>
> I get numbers that are not even closed to the truth!<br>
> Thank you!<br>
          <br>
Nena, in suggesting the "backwards" technique on the integral pseries 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 <br>
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. <br>
<br>
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 nontrivial task, not to be attained in such an easy manner.<br>
<br>
Roger Stafford

Tue, 01 Mar 2011 02:26:05 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822315
Nena Kabranski
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?<br>
Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.<br>
Could you help me in this analysis?

Tue, 01 Mar 2011 04:43:05 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822325
Roger Stafford
"Nena Kabranski" wrote in message <ikhlft$3h1$1@fred.mathworks.com>...<br>
> 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?<br>
> Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.<br>
> Could you help me in this analysis?<br>
          <br>
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 53bit 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.<br>
<br>
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 <br>
<br>
bessef(3)=2*1*bessef(2)bessef(1)<br>
<br>
bessef(3) will accordingly be too high by .0000000002. Then pursuing this further we get<br>
<br>
bessef( 4) too high by .0000000007<br>
bessef( 5) too high by .0000000040<br>
bessef( 6) too high by .0000000313<br>
bessef( 7) too high by .0000003090<br>
......<br>
bessef(12) too high by .2915473799<br>
<br>
In this rapid buildup 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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
Roger Stafford

Tue, 01 Mar 2011 14:02:05 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822424
Matt J
"Roger Stafford" wrote in message <ikhtgp$7e0$1@fred.mathworks.com>...<br>
><br>
> 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.<br>
===============<br>
<br>
@Nena, <br>
<br>
...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.

Tue, 01 Mar 2011 15:35:28 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822448
Steven_Lord
<br>
<br>
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in <br>
message news:ikhtgp$7e0$1@fred.mathworks.com...<br>
> "Nena Kabranski" wrote in message <ikhlft$3h1$1@fred.mathworks.com>...<br>
>> I guess I am not trying to solve those problems… I am trying to <br>
>> understand in which cases the MatLab is making mistakes. I cannot see a <br>
>> problem here, because the numbers are close to each other and they are <br>
>> not extremely small. Or are they?<br>
>> Abs (J(n))<=1, because it is an area of a function that is at most 1/pi.<br>
>> Could you help me in this analysis?<br>
>           <br>
> Nena, I am trying to get you to see that it isn't simply Matlab that is <br>
> making the errors you have found. With the algorithm you have used, any <br>
> computing system possessing the same 53bit accuracy limitation would make <br>
> the same kinds of errors. This is inherent in the very nature of digital <br>
> computing. It is the algorithm being used that must be considered to be <br>
> inappropriate.<br>
<br>
*snip Roger's excellent example*<br>
<br>
In fact, there's an entire branch of mathematics dedicated to this type of <br>
problem.<br>
<br>
<a href="http://en.wikipedia.org/wiki/Numerical_analysis">http://en.wikipedia.org/wiki/Numerical_analysis</a><br>
<br>
You might also find this Cleve's Corner article interesting and informative:<br>
<br>
<a href="http://www.mathworks.com/company/newsletters/news_notes/clevescorner/winter02_cleve.html">http://www.mathworks.com/company/newsletters/news_notes/clevescorner/winter02_cleve.html</a><br>
<br>
as well as this oftencited paper:<br>
<br>
<a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.102.244">http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.102.244</a><br>
<br>
 <br>
Steve Lord<br>
slord@mathworks.com<br>
To contact Technical Support use the Contact Us link on <br>
<a href="http://www.mathworks.com">http://www.mathworks.com</a>

Wed, 02 Mar 2011 03:01:23 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#822580
Nena Kabranski
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. :)<br>
Sincerely,<br>
Nena<br>
<br>
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.

Wed, 23 Mar 2011 22:53:02 +0000
Re: where is the mistake
http://www.mathworks.com/matlabcentral/newsreader/view_thread/303553#827008
Mark Shore
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ikdrci$jch$1@fred.mathworks.com>...<br>
<br>
> <br>
> I could not resist trying out my new toolbox in MATLAB.<br>
> <br>
> FPF is a class that allows you to work in floating point<br>
> arithmetic with a user specified precision. So here, with<br>
> 50 digits of precision in all computations...<br>
> <br>
> intxex=fpf(1,50);<br>
> y=1;<br>
> n = 20;<br>
> e1 = fpf('e',50);<br>
> for i=2:n<br>
> intxex=e1i*intxex;<br>
> y=intxex<br>
> end<br>
> <br>
<br>
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.