Getting different results from function handle & syms for a same equation. How to avoid it?
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
Matt J
2021-7-30
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 个评论
Matt J
2021-7-30
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
2021-7-30
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 个评论
Walter Roberson
2021-8-2
编辑:Walter Roberson
2021-8-2
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
2021-8-2
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!