Unrealistic results from integration

1 次查看（过去 30 天）
Svante Monie
Svante Monie 2015-7-23
编辑： Walter Roberson 2015-7-27
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

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2015-7-23
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 个评论
Svante Monie
Svante Monie 2015-7-27
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

请先登录，再进行评论。

Titus Edelhofer
Titus Edelhofer 2015-7-23
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 个评论
Svante Monie
Svante Monie 2015-7-27
Hello Titus
Thanks for the cleaning up tips!
Looks much nicer, but unfortunately the issue persists.
//Svante

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by