Can I solve the integral given below using Matlab?

$$\frac{C(J)}{C_0} = \frac{2e^J}{\pi}\int\limits_{0}^{\infty} \frac{e^{\frac{\gamma}{2}\left[1 - \sqrt{\rho}\cos(\theta/2)\right]}}{a_1^2 + a_2^2} [a_1\cos(ZJ-w) + a_2\sin(ZJ - w)]\, dZ$$ where \begin{align} &\theta = \arctan(v/u)\\ &u = 1 + \frac{4}{\gamma}\left[1 + \frac{ba + a(1 + Z^2)}{(1 + b)^2 + Z^2}\right]\\ &v = \frac{4Z}{\gamma}\left[1 + \frac{ab}{(1 + b)^2 + Z^2}\right]\\ &\rho = \sqrt{u^2 + v^2}\\ &b = \frac{af}{1 - f}\\ &a_1 = 1 + \sqrt{\rho}\cos(\theta/2) - Z\sqrt{\rho}\sin(\theta/2)\\ &a_2 = Z\left[1 + \sqrt{\rho}\cos(\theta/2)\right] + \sqrt{\rho}\sin(\theta/2)\\ &w = \frac{\gamma}{2}\sqrt{\rho}\sin(\theta/2) \end{align}

Relation between $Z$ and $J$ is as follows: $$Z = \frac{af}{1 - f} (J-y) \enspace .$$ $y$ varies in $[0, 1]$; $C_0=1$; $a= 0.065$; $f=0.7$. I am new to Matlab. Any help will be highly appreciated. I ran the following code (gamma is not y and gamma=m=60, J=1/0.7, O=theta, C=C(J)/C0), but matlab is busy from last 12-13 hours!!

syms z
a = 0.065;
f = 0.7;
b = a*f/(1-f);
u = 1+((4/m)*(1+(b*a+a*(1+z.^2))/((1+b).^2)+z.^2));
v = (4*z/m)*(1+(a*b/(1+b).^2+z.^2));
p = sqrt(u.^2+v.^2);
O = atan(v/u);
a1 = 1+sqrt(p)*cos(O/2)-z*sqrt(p)*sin(O/2);
a2 = z*(1+sqrt(p)*cos(O/2))+sqrt(p)*sin(O/2);
w = (m/2)*sqrt(p)*sin(O/2);
fun = (exp(m/2*(1-sqrt(p)*cos(O/2)))/(a1.^2+a2.^2))*(a1*cos((z/0.7)-w)+a2*sin((z/0.7)-w));
I = int(fun, z, 0, inf);
C = 2*exp(1/0.7)*I/pi;
if C>1
  • 1
    $\begingroup$ Have you tried the integral function, or one of the other integration functions in matlab, and did it work? If it didn't work, can you say how it failed? $\endgroup$ – Kirill May 15 '15 at 19:19
  • $\begingroup$ I just edited your equations, but it seems that what you call $y$ is $\gamma$. Am I correct? $\endgroup$ – nicoguaro May 15 '15 at 19:26
  • 1
    $\begingroup$ @user16362 please confirm if the edited equations are correct $\endgroup$ – Stefano M May 15 '15 at 20:57
  • $\begingroup$ @nicoguaro, sorry I forgot to mention y is not γ. γ=60. I tried to run the code mentioned in answer with J=1/0.7, but matlab is busy from last 12 hours!! $\endgroup$ – user16362 May 16 '15 at 10:58
  • 1
    $\begingroup$ If you want to solve this numerically, then you probably should be using symbolic math and the int function. Try integral as suggested by @Kiril: integral(matlabFunction(simplify(fun)),0,Inf). Although with an oscillatory integrand, you may want to try alternative quadrature schemes to check your answer (the exponential term may provide sufficient damping at large values however). See this question $\endgroup$ – horchler May 16 '15 at 20:37

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.