Unrealistic results from integration

2 views (last 30 days)
Svante Monie
Svante Monie on 23 Jul 2015
Edited: Walter Roberson on 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

Answers (2)

Torsten
Torsten on 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 Comment
Svante Monie
Svante Monie on 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

Sign in to comment.


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

Sign in to comment.

Categories

Find more on Characters and Strings in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!