How do I integrate an erfc (complementary error function)?

6 次查看(过去 30 天)
I'm trying to calculate following in Matlab:
U = 1/pi * integral from [0, 1/Tao] {erfc(gamma*f(x)/sqrt(Tao - f(x))} dx where f(x) = 3 * sin(x) - x * cos(x)/(sin(x))^3 gamma, Tao = variables
I've tried many different ways, but none works. The last attempt looks like this:
gamma = 0.7993;
TAO = [0,0.6127,1.2255,1.8382,2.4509,3.0636,3.6764,4.2891,4.9018,5.5145,6.1273];
f0TAO = zeros();
for i=2:length(TAO) f0TAO(i) = (3*((sin(TAO(i))-(TAO(i)*cos(TAO(i))))/(sin(TAO(i)))^3));
end
f0Tao = 1./f0TAO;
Q = zeros();
for i=2:length(Tao);
c = TAO(i);
f = erfc(gamma.*((3.*((sin(x)-(x.*cos(x)))./(sin(x))^3))./sqrt(c-(3.*((sin(x)-(x.*cos(x)))./(sin(x))^3)))));
Q(i) = 1/pi*int(f,0,f0Tao(i));
end
This resulted in an error like this:
The following error occurred converting from sym to double: Error using mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.
If the input expression contains a symbolic variable, use the VPA function instead.
What should I do? Can anyone give me some brilliant clues?
Thanks in advance!

采纳的回答

Torsten
Torsten 2015-7-15
Q(i)=1/pi*vpa(int(f,0,f0Tao(i)),5);
But I guess the integration will be difficult because you divide by sin(x)^3 which is zero at x=0.
Additionally, the sqrt will become complex-valued.
Best wishes
Torsten.
  8 个评论
Torsten
Torsten 2015-7-15
If
f=@(x)erfc(gamma.*((3.*((sin(x)-(x.*cos(x)))./(sin(x))^3))./sqrt(c-(3.*((sin(x)-(x.*cos(x)))./(sin(x))^3)))));
f(2);
still gives an error, check the value of the argument to the erfc function and whether gamma and c are real and not symbolic.
Best wishes
Torsten.
Svante Monie
Svante Monie 2015-7-15
编辑:Svante Monie 2015-7-15
Now I have found the cause of the error. In the denominator I have sqrt(c-(3.*((sin..., this gives complex values and thats what MatLab don't like. It needs _real_ input, as stated. When I put abs() around the expression under the sqrt, erfc delievers! Question is then of course, if the integral is valid for my intentions...

请先登录,再进行评论。

更多回答(0 个)

产品

Community Treasure Hunt

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

Start Hunting!

Translated by