Unrealistic results from integration

1 Ansicht (letzte 30 Tage)
Svante Monie
Svante Monie am 23 Jul. 2015
Bearbeitet: Walter Roberson am 27 Jul. 2015
I'm trying to implement some mathematic formulas to MatLab, but when I excecute the program with parameters accordning to a referens case I get very unlikely result. Somewhere something takes the wrong turn. Basically I'm trying to perform this integral:
U_out(Tao) = 1/pi * integral (from 0 to 1/f_0(Tao)) {erfc(gamma * (f_0(s)/sqrt(Tao-f_0(s))))} ds
where f_0(s) = 3 * ((sin(s) -s*cos(s))/sin^3(s))
Tao is a vector [0 1.51 3.02 4.53 6.04 7.55]
gamma = 0.457
s = 197.23
1/f_0(Tao) = 1*10^-3 * [inf -0.0041 -0,0164 -0.0370 -0.0656 -0.1029]
According to the reference case I'm supposed to get following U_out's:
@ Tao = 2.268 => U_out appr 0.18
@ Tao = 4.532 => U_out appr 0.29
@ Tao = 7.553 => U_out appr 0.37
The results I get is unfortunately negative and four potenses too small Also I tend to get the same result regardless if I excecute the integral with the variable 's' or with 'x'.
Some exctract from the code:
gamma = 0.456915136131553
Tao = 0 1.510551431477735 3.021102862955471 4.531654294433206 6.042205725910941 7.552757157388677
f0Tao = 1.0e-03 * Inf -0.004115574605475 -0.016462509923427 -0.037041490055367 -0.065853821127324 -0.102901682460165
f0s = -2.169971152489210e+04
u = 2 % Tao must be larger then 1
U_out = zeros();
for i=u:length(Tao);
c = Tao(i);
f = @(x) erfc(gamma.*((3.*((sin(x)-(x.*cos(x)))./(sin(x)).^3))./sqrt(c-(3*((sin(x)-(x.*cos(x)))./(sin(x)).^3)))));
U_out(i) = 1/pi*integral(f,0,f0Tao(i));
%f = @(x) erfc(gamma.*(f0s./sqrt(c-f0s)));
%U_out(i) = 1/pi*integral(f,0,f0Tao(i),'ArrayValued',true);
end

Antworten (2)

Torsten
Torsten am 23 Jul. 2015
You should check
1/f_0(Tao) = 1*10^-3 * [inf -0.0041 -0,0164 -0.0370 -0.0656 -0.1029]
again.
I get for the first two values
[1 0.3657 ...]
Best wishes
Torsten.
  1 Kommentar
Svante Monie
Svante Monie am 27 Jul. 2015
Hello Torsten, sorry for my late respons. I guess you used the radius function of sine and cosine for 1/f_0(Tao)?
When doing that I get similar values for 1/f_0(Tao): [NaN 0.3654 0.000186 1.876...] But these don't follow any logical trend (it's all about temperature decline so it should be some sort of logical trend i.e. linear or exponentiell).
Furthermore this gives complex values for the function erfc as I get negative values in the sqrt and erfc can't handle that...
If I use absolute value I get pass that problem, but the values I then get is completely rubbish:
U_out [0 0.0405 0.00004 0.29 0.00018 -0.057]
The only way to get erfc to accept the input is by using the degree function for sine and cosine, but then I'm back at square one...
//Svante

Melden Sie sich an, um zu kommentieren.


Titus Edelhofer
Titus Edelhofer am 23 Jul. 2015
Hi,
I'm not getting the values that are supposed to be the correct values, but here is a slightly cleaner version of your code as starting point:
% define f_0 as a function:
f0 = @(s) 3 .* (sin(s) -s.*cos(s))./(sin(s).^3);
% do computation for one Tao and gamma:
Tao = 2.268;
gamma = 0.456915136131553;
% define the function f using f_0 from above
f = @(s) erfc(gamma * f0(s)./sqrt(Tao-f0(s)));
% compute the integral (note, the right boundary should be 1/f0(Tao)?
U_0 = 1/pi * integral(f, 0, 1/f0(Tao))
Hope this helps,
Titus
  1 Kommentar
Svante Monie
Svante Monie am 27 Jul. 2015
Hello Titus
Thanks for the cleaning up tips!
Looks much nicer, but unfortunately the issue persists.
//Svante

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Fourier Analysis and Filtering finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by