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

5 visualizzazioni (ultimi 30 giorni)
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!

Risposta accettata

Torsten
Torsten il 15 Lug 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 Commenti
Torsten
Torsten il 15 Lug 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 il 15 Lug 2015
Modificato: Svante Monie il 15 Lug 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...

Accedi per commentare.

Più risposte (0)

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by