Getting different results from function handle & syms for a same equation. How to avoid it?
3 views (last 30 days)
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.
0 Comments
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
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
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.
Walter Roberson
on 2 Aug 2021
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!