Major confusion for loops - precision not acquired

3 次查看（过去 30 天）
DM
DM 2015-2-6
评论： Star Strider 2015-2-6
I have a loop that sums the values of function called pdfomegat
omegat= -1.9999:0.01:+1.9999;
sum1=0;
for j=1:length(omegat)
func(j)=pdfomegat(j)*0.01;
sum1=sum1+func(j);
end
c1=sum1;
My second loop is slightly different loop and is as follows
sum2=0;
syms t
for j=1:length(omegat)
func2(j)=t*pdfomegat(j)*0.01;
sum2=sum2+func(j);
end
c2=sum2;
I know for fact that c1 is 1 .. However notice that in the second loop I change nothing but multiply by symbolic character t. Howver I get that
c2=0.99999464659732599597363744692302*t.
Why isn't c2=t?
Is there anyway I can fix it, I can't use the round function because it symbolic.. Any thoughts?

请先登录，再回答此问题。

回答（3 个）

A Jenkins
A Jenkins 2015-2-6
Play around with vpa() and digits() and make sure to read the help files on Precision Control.
  2 个评论
DM
DM 2015-2-6
Is there any way I can extract the number from the symbolic expression and round it off? My problem requires high accuracy.
A Jenkins
A Jenkins 2015-2-6
>> syms t
>> c2=0.99999464659732599597363744692302*t;
>> cf=vpa(c2,8)
cf =
0.99999465*t
>> cf=vpa(c2,4)
cf =
1.0*t

请先登录，再进行评论。

John D'Errico
John D'Errico 2015-2-6
编辑：John D'Errico 2015-2-6
No. You cannot make a double precision number more precise than it is. There are simple many numbers that are not exactly representable using a double. For example, 0.1 is not so, since doubles are stored in binary form.
Anyway, my very good bet is that c1 is NOT 1, but that it was displayed in the command window as 1, using format short. So I think what you know for a fact is actually completely false. As well, it is entirely possible that you were using a single precision format to store pdfomegat, since c1 was so inaccurately NOT equal to 1.
However, you can do many things to improve your work. For example, in this line, rather than multiplying by 0.01, divide by 100, an integer. You do this because 0.01 is not stored exactly as a double, but 100 is.
As for the code itself, write it more efficiently. Learn to use vectors, and simple tools like sum.
Your first loop does nothing more than multiply the vector pdfomegat by 0.01, and then sum it.
omegat= -1.9999:0.01:+1.9999;
func1 = pdfomegat/100;
c1 = sum(func1);
You can also totally replace that second loop simply as
syms t
func2 = t*func1;
c2 = t*c1;
As for making the sums themselves more accurate, you could convert the vector pdfomegat into a symbolic form (or my own HPF) but since they were originally stored as doubles, you simply won't gain much at all there, unless you create those numbers in a form that allows you to store them more accurately in the first place.
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2015-2-6
But you are guessing completely WRONG.
That you get 1.0000 on the command line proves my point completely! It was displayed under format short rules.
help format
See what
format long g
would do instead. Then you would learn that your number was not in fact exactly 1, but slightly off. When the symbolic toolbox gets that number, it displays the number more accurately.
As for the code I showed, I would still contend that you need to learn to use vectors, that you could do what you are doing more efficiently and more cleanly. I cannot know that positively, but given the code you did show, it seems a very probable conclusion.
John D'Errico
John D'Errico 2015-2-6
Any rounding that you do AFTER the fact merely hides the numerical issues. It solves no problem. Instead, you need to show more clearly what you are doing, and do that computation more accurately.

请先登录，再进行评论。

Star Strider
Star Strider 2015-2-6
It’s not straightforward, but you can round the coefficient:
syms t
c2=0.99999464659732599597363744692302*t; % Original Expression
cf = double(vpa(coeffs(c2))) % Get Coefficient
cf = fix(cf * 1E+8)/1E+8 % Round To 8 Digits
  2 个评论
DM
DM 2015-2-6
Thanks, but what if the function was not polynomial, my function is of the folllowing form
func2= 1/(t^-3+1).*func1*g;
I am sorry I really tried to keep everything very simple. Coeffs doesnt work on this expression.
Star Strider
Star Strider 2015-2-6
I would suggest that ‘keeping everything simple’ is not appropriate. We Answer the Questions posted, so if you do not state your actual problem, you will likely not get an Answer that will solve it.
I have no idea what ‘func1’ and ‘g’ are, so I can go no further. You might use the numden function to separate the numerator and denominator, then use coeffs or whatever function is appropriate. Consider also matlabFunction or similar functions if you want to use your function in other applications.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Conversion Between Symbolic and Numeric 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by