Getting different results from function handle & syms for a same equation. How to avoid it?

1 visualización (últimos 30 días)
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)),
r = 5×1
0 7.9332 145.4839 217.7066 217.7066
>> 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),
R = 5×1
0 7.9332 156.3824 197.9575 249.7273
The resulst from syms looks correct.
kindly give me a suggestion to avoid the resuls deviation in function handle.

Respuestas (2)

Matt J
Matt J el 30 de Jul. de 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 comentarios
Matt J
Matt J el 30 de Jul. de 2021
Editada: Matt J el 30 de Jul. de 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
Matt J el 30 de Jul. de 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

Iniciar sesión para comentar.


Walter Roberson
Walter Roberson el 30 de Jul. de 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 comentarios
Walter Roberson
Walter Roberson el 2 de Ag. de 2021
Editada: Walter Roberson el 2 de Ag. de 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
Walter Roberson el 2 de Ag. de 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)
P = 
[S,T] = coeffs(P, u, 'all')
S = 
T = 
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])
SF = function_handle with value:
@(c1_12,c2_21,c3_12,c4_21,c5_12,c6_21,c7_12,c8_21,c9_12)[-c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21.*c7_12.*c8_21.*c9_12,c9_12.*(c8_21.*(c7_12.*(c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21.*c7_12,-c9_12.*(c8_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c7_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)-c7_12.*(c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)-c1_12.*c2_21.*c3_12.*c4_21.*c5_12,c9_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c8_21.*(c1_12+c3_12+c5_12+c7_12)+c4_21.*(c1_12+c3_12))+c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c7_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12,-c1_12-c3_12-c5_12-c7_12-c9_12,0.0]
func2str(SF)
ans = '@(c1_12,c2_21,c3_12,c4_21,c5_12,c6_21,c7_12,c8_21,c9_12)[-c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21.*c7_12.*c8_21.*c9_12,c9_12.*(c8_21.*(c7_12.*(c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21)+c1_12.*c2_21.*c3_12.*c4_21.*c5_12.*c6_21.*c7_12,-c9_12.*(c8_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c7_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)-c7_12.*(c6_21.*(c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12)+c1_12.*c2_21.*c3_12.*c4_21)-c1_12.*c2_21.*c3_12.*c4_21.*c5_12,c9_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c8_21.*(c1_12+c3_12+c5_12+c7_12)+c4_21.*(c1_12+c3_12))+c5_12.*(c1_12.*c2_21+c4_21.*(c1_12+c3_12))+c7_12.*(c1_12.*c2_21+c6_21.*(c1_12+c3_12+c5_12)+c4_21.*(c1_12+c3_12))+c1_12.*c2_21.*c3_12,-c1_12-c3_12-c5_12-c7_12-c9_12,0.0]'
Basically... avoid the inaccuracy problems by skipping the fitting since it is not needed.

Iniciar sesión para comentar.

Community Treasure Hunt

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

Start Hunting!

Translated by