0
$\begingroup$

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
m=60;
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;
display(C)
plot(z,C)
xlabel('Z')
ylabel('concentration')
if C>1
    stop
end
$\endgroup$
  • 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$ @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
  • $\begingroup$ Do you want to integrate that equation analytically? MATLAB is not the best tool for that aim. $\endgroup$ – nicoguaro May 16 '15 at 12:11
  • $\begingroup$ @nicoguaro I want to solve it numerically! should I change the upper limit of integral to some finite number? $\endgroup$ – user16362 May 16 '15 at 18:22
  • 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.