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

6 views (last 30 days)
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!

Accepted Answer

Torsten
Torsten on 15 Jul 2015
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 Comments
Torsten
Torsten on 15 Jul 2015
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 on 15 Jul 2015
Edited: Svante Monie on 15 Jul 2015
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...

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!