Unrealistic results from integration

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

Answer 1

[Torsten]

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.

[Svante Monie]

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

Answer 2

[Titus Edelhofer]

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

[Svante Monie]

Hello Titus

Thanks for the cleaning up tips!

Looks much nicer, but unfortunately the issue persists.

//Svante