Getting different results from function handle & syms for a same equation. How to avoid it?
Show older comments
I am in the process of replacing syms with the function handle.
I have created solution approaches for my application which is finding roots of a polynomial eqation.
I am getting different results from syms & function handle.
The results deviation increases with increase in the polynomial order of the equation.
I am sharing the both the scripts for the reference.
>> Function handle
f = @(u) [1 -u*0.112;0 1]*[1 0;0.000945465 1]*[1 -u*0.0358;0 1]*[1 0;0.106859097 1]*[1 -u*0.057780033;0 1]*[1 0;0.000473198 1]*[1 -u*0.243621246;0 1]*[1 0;9.98448E-05 1]*[1 -u*33.60710411;0 1];
%finding the roots of the eqation at final matrix element (1,2)%
for u = 0:5
P=f(u);
V(u+1,1)=P(1,2);
M(u+1,:) = [u^5,u^4,u^3,u^2,u,1];
end
coeffs = M\V; % the co-efficients of the polynomial at matrix element (1,2)
r = roots(coeffs);
r=abs(r);
r=sort(sqrt(r)),
>> With syms
syms u
f = [1 -u*0.112;0 1]*[1 0;0.000945465 1]*[1 -u*0.0358;0 1]*[1 0;0.106859097 1]*[1 -u*0.057780033;0 1]*[1 0;0.000473198 1]*[1 -u*0.243621246;0 1]*[1 0;9.98448E-05 1]*[1 -u*33.60710411;0 1];
%finding the roots of the eqation at final matrix element (1,2)%
P=f(1,2);
R=vpasolve(P);
R=double(R);
R=sqrt(R),
The resulst from syms looks correct.
kindly give me a suggestion to avoid the resuls deviation in function handle.
Answers (2)
Matt J
on 30 Jul 2021
The results are different because the polynomial coefficients are different:
>> sym2poly(P)-coeffs'
ans =
1.0e-04 *
0.0000 0.0000 -0.0000 -0.3802 0.0539 0
3 Comments
Vinothkumar Sethurasu
on 30 Jul 2021
No, in addition to the things that Walter said, the condition number of your Vandermonde matrix is quite high.
>> cond(M)
ans =
5.7689e+04
This will make the results of M\V more error prone.
Matt J
on 30 Jul 2021
Notice also that if you normalize the leading coefficient to 1, you can see that your coefficient values nearly push the limit of double float precision:
>> coeffs/coeffs(1)
ans =
1.0e+15 *
0.0000
-0.0000
0.0000
-0.0467
2.9237
0
Walter Roberson
on 30 Jul 2021
Remember that when you mix floating point numbers into a symbolic expression, that MATLAB uses the default conversion of floating point to symbolic numbers, which is the 'r' (rational) conversion.
sr = sym(0.0358)
sd = sym(0.0358,'f')
sr - sd
sr is what will be converted into by default, 179/5000 exactly, which is 358/10000 .
sd is the exact binary fraction that is stored for double precision 0.0358 . The denominator is 2^54
The two differ...
These sorts of little differences accumulate quickly.
I am getting different results from syms & function handle.
The results deviation increases with increase in the polynomial order of the equation.
-u*33.60710411
What kind of experiment are you doing, such that you are able to measure one of the values to 10 digits of precision, but you are only able to measure 0.0358 to 3 digits of precision? And why are you expecting more than a small number of digits of precision on the roots when you have inputs that are only 3 digits of precision?
Remember that as far as Science is concerned, when you wrote 0.0358 then you mean "a number whose exact value is not known, but which is between 3575/100000 (inclusive) and 3585/100000 (exclusive). With such a wide range, how can you expect to get "exact" solutions to the roots?
8 Comments
Vinothkumar Sethurasu
on 30 Jul 2021
Edited: Vinothkumar Sethurasu
on 30 Jul 2021
Vinothkumar Sethurasu
on 2 Aug 2021
Edited: Vinothkumar Sethurasu
on 2 Aug 2021
Walter Roberson
on 2 Aug 2021
My suggestion to avoid the result deviation is
- use the maximum known precision for each number. Include all trailing zeros: they might not make a difference in the calculation, but they make a difference in documentation that you know what you are doing
- If the numbers in your equations are derived numbers, such as if you have replaced 1/1057.68 with 0.000945465 then put in the original numbers rather than the derived numbers, so that you get the highest precision
- If you can, use the Symbolic toolbox.
- When you use the Symbolic Toolbox, rewrite every floating point number as a rational number. For example rewrite 0.106859097 as sym(106859097)/10^10. Take advantage of any information you know about the numbers to write the best possible rational version
- When you use the Symbolic Toolbox, rewrite every number involving an exponent as a product of a rational and a power of 10. For example, 9.98448E-05 would be written as sym(998448)*10^(-10)
- When you use the Symbolic Toolbox, simplify() the expressions
- If you must use floating point calculations, then if possible write the expressions in terms of the Symbolic Toolbox and then use matlabFunction()
If you do not have the symbolic toolbox, then all you can do is use maximum precision in how you write your numbers, together with any hand calculation to re-arrange the terms to increase precision.
There is no way to tell MATLAB to use "quad" precision floating point numbers.
for u = 0:5
P=f(u);
V(u+1,1)=P(1,2);
M(u+1,:) = [u^5,u^4,u^3,u^2,u,1];
end
I suggest you extract the P(1,2) for each u value, but after that use polyfit() instead of building the Vandermode matrix yourself.
Vinothkumar Sethurasu
on 2 Aug 2021
format long g
Q = @(v) sym(v);
f = @(u) [1 -u*Q(112)/10^3;0 1]*[1 0;Q(945465)/10^9 1]*[1 -u*Q(358)/10^4;0 1]*[1 0;Q(106859097)/10^9 1]*[1 -u*Q(57780033)/10^9;0 1]*[1 0;Q(473198)/10^9 1]*[1 -u*Q(243621246)/10^9;0 1]*[1 0;Q(998448)/10^10 1]*[1 -u*Q(3360710411)/10^8;0 1];
syms U
p = f(U)
F = simplify(p(1,2))
V = subs(F,U,0:5)
polyfit(0:5,double(V), 5)
accurate_values = sym2poly(F)
Notice that F is already in the form of a polynomial of the appropriate degree, so all that really needs to be done is sym2poly() to get the accurate coefficients.
Now let's go back and see what we can get from the original function handle:
ff = @(u) [1 -u*0.112;0 1]*[1 0;0.000945465 1]*[1 -u*0.0358;0 1]*[1 0;0.106859097 1]*[1 -u*0.057780033;0 1]*[1 0;0.000473198 1]*[1 -u*0.243621246;0 1]*[1 0;9.98448E-05 1]*[1 -u*33.60710411;0 1];
syms U
pf = ff(U);
Ff = pf(1,2);
less_accurate_values = sym2poly(Ff)
accurate_values - less_accurate_values
max(abs(ans))
Those are pretty accurate actually.
So... if you work symbolically then you do not need to do any fitting at all.
Vinothkumar Sethurasu
on 2 Aug 2021
Walter Roberson
on 2 Aug 2021
Edited: Walter Roberson
on 2 Aug 2021
You cannot do that -- not without expanding out your values algebraically. Matt pointed out that you are pushing the limits of double precision.
syms u c1_12 c2_21 c3_12 c4_21 c5_12 c6_21 c7_12 c8_21 c9_12 real
ff = [1 -u*c1_12;0 1]*[1 0;c2_21 1]*[1 -u*c3_12;0 1]*[1 0;c4_21 1]*[1 -u*c5_12;0 1]*[1 0;c6_21 1]*[1 -u*c7_12;0 1]*[1 0;c8_21 1]*[1 -u*c9_12;0 1];
P = ff(1,2)
[S,T] = coeffs(P, u, 'all')
So your matrix of coefficients that you are looking for is in S, and it is in the order of u^5 down to constant. Instead of doing the calculation that you are doing, you can construct S as a function that takes in the numeric constants and returns the coefficients.
SF = matlabFunction(S, 'vars', [c1_12 c2_21 c3_12 c4_21 c5_12 c6_21 c7_12 c8_21 c9_12])
func2str(SF)
Basically... avoid the inaccuracy problems by skipping the fitting since it is not needed.
Categories
Find more on Number Theory in Help Center and File Exchange
on 30 Jul 2021
on 2 Aug 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
