Problem with plotting an anonymous function

calka1 = @(X)(integral(@(x)((x.*(-1.15517395637196e-2)+x.^2.*5.25774286994568e-4-x.^3.*1.089934292266744e-5-x.^4.*1.230239476736945e-7+x.^5.*6.42213403953954e-10+x.^6.*1.550438107746178e-11+x.^7.*3.670339539265712e-15-x.^8.*1.053763407268398e-15-x.^9.*1.00711819287645e-18+x.^10.*3.225336064282519e-20+x.^11.*1.85216906506446e-23-x.^12.*3.481067131642152e-25+2.92686659979064e-2)./(X-x)), -180, 180))
When I try to plot an anonymous function with
fplot(calka1, [-180 180])
I get
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is
4.1e+00. The integral may not exist, or it may be difficult to approximate numerically to the requested
accuracy.
> In integralCalc/iterateScalarValued (line 372)
In integralCalc/vadapt (line 132)
In integralCalc (line 75)
In integral (line 88)
In @(X)(integral(@(x)((x.*(-1.15517395637196e-2)+x.^2.*5.25774286994568e-4-x.^3.*1.089934292266744e-5-x.^4.*1.230239476736945e-7+x.^5.*6.42213403953954e-10+x.^6.*1.550438107746178e-11+x.^7.*3.670339539265712e-15-x.^8.*1.053763407268398e-15-x.^9.*1.00711819287645e-18+x.^10.*3.225336064282519e-20+x.^11.*1.85216906506446e-23-x.^12.*3.481067131642152e-25+2.92686659979064e-2)./(X-x)),-180,180))
In fplot (line 128)
spammed into my face and the plot I recieve is jumping from 300 to -500 like crazy.
Am I doing something wrong?

댓글 수: 3

Stephen23
Stephen23 2020년 8월 18일
편집: Stephen23 2020년 8월 18일
"Am I doing something wrong?"
Numeric integration of a degree 12 polynomial... aka "noise amplification". I would be surprised if an untailored numeric integration method could provide a stable result.
But it has no problem with integrating this one
int(7.30561394644000E+02 + 2.92686659979064E-02*x -5.77586978185980E-03*x^2 + ...
1.75258095664856E-04*x^3 -2.72483573066686E-06*x^4 -2.46047895347389E-08*x^5 + ...
1.07035567325659E-10*x^6 + 2.21491158249454E-12*x^7 + 4.58792442408214E-16*x^8 ...
-1.17084823029822E-16*x^9 -1.00711819287645E-19*x^10 + 2.93212369480229E-21*x^11 + ...
1.54347422088705E-24*x^12 -2.67774394741704E-26*x^13, -180, 180)
as it gives me a good result. These functions are taken from CurveExpert which made it out of ship hull's shape.
The function in the main post is actually a derivative from this one and divided by (X-x). Is it possible that it's too much? I mean I don't really have a choice with that thing.
Thank you for the comment tho.
Stephen23
Stephen23 2020년 8월 18일
편집: Stephen23 2020년 8월 18일
"...as it gives me a good result."
Perhaps.
Or then again, maybe not.
Given the small** coefficient values (down to 1e-26) paired with large** values calculated from x (up to 1e30) your numeric data cover quite a wide range, and operations on them could easily result in some loss of information or some other similar floating point effects. Certainly without checking the numeric validity of these operations (from a binary floating number point of view) caution in trusting the output would be wise.
** relative to each other, in a binary floating number sense.

댓글을 달려면 로그인하십시오.

 채택된 답변

Hubert P
Hubert P 2020년 8월 20일

0 개 추천

The reason is the singularity (X-x) in the denominator.

추가 답변 (1개)

Rik
Rik 2020년 8월 18일
편집: Rik 2020년 8월 18일

1 개 추천

Let's first make the terms more readable:
calka1 = @(X)(integral(@(x)((...
x.*(-1.15517395637196e-2)...
+x.^2.*5.25774286994568e-4...
-x.^3.*1.089934292266744e-5...
-x.^4.*1.230239476736945e-7...
+x.^5.*6.42213403953954e-10...
+x.^6.*1.550438107746178e-11...
+x.^7.*3.670339539265712e-15...
-x.^8.*1.053763407268398e-15...
-x.^9.*1.00711819287645e-18...
+x.^10.*3.225336064282519e-20...
+x.^11.*1.85216906506446e-23...
-x.^12.*3.481067131642152e-25...
+2.92686659979064e-2...
)./(X-x)), -180, 180, 'ArrayValued',1));
calka1(0)
When we run this we already get a warning that the bound on the error is 1e-1. That seems quite large. Can you see an indication of why this happens?
Welcome to the wondrous world of floats. You have very large exponents which are paired with tiny factors.
180^12 %about 1e27
flintmax %about 9e15
You can't even store every integer at your outer bounds.
My rule of thumb: always make sure your exponents don't span more than about 10, but at most eps:
eps('double') % 2.2204e-16

댓글 수: 1

I used the function
@(X)(integral(@(x)((sin(x.*1.006116719101209e-2-1.066548307975591e-1).*(-6.137229496083432))./(X-x)),-180,180,'ArrayValued',1))
And the results are quite the same + the error is the same (or even bigger).
Any ideas?

댓글을 달려면 로그인하십시오.

카테고리

제품

릴리스

R2015a

질문:

2020년 8월 18일

답변:

2020년 8월 20일

Community Treasure Hunt

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

Start Hunting!

Translated by