Unrealistic results from integration
1 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(2 个)
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.
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Characters and Strings 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!