How to integrate such complicated expression in matlab?

1 view (last 30 days)
P =
where is Gamma function and the next Gamma is upper incomplete Gamma and are all constants.
My query is how to integrate this in Matlab.
I will appreciate any help regarding this.

Accepted Answer

Walter Roberson
Walter Roberson on 28 Sep 2021
syms A B gamma_s k m sigma_e sigma_r y real
inner = (gamma(1+m) + igamma( 1+m, (gamma_s*(1+y)-1) / (A*sigma_r^2) )) .* y.^k .* exp(-y/(B*sigma_e)^2)
inner = 
P = int(inner, y, 0, inf)
P = 
You are not going to get a closed form solution for symbolic k (and probably not for symbolic m either.)
  2 Comments
Walter Roberson
Walter Roberson on 28 Sep 2021
If you subs() in numeric k and m, then you will get the result as a limit() as y approaches infinity. If you add the assumption that the constants are all positive, then the limit will be resolved.
For example with k = 3, m = 4 and simplify() and collect with regards to sigma_e and gamma_s then
(144*B^24*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^24*gamma_s^8 + 1152*A*B^22*sigma_r^2*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^22*gamma_s^7 + 4032*A^2*B^20*sigma_r^4*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^20*gamma_s^6 + 8064*A^3*B^18*sigma_r^6*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^18*gamma_s^5 + 6*B^16*sigma_e^16*gamma_s^8 + (120*A*B^16*sigma_r^2 - 24*B^16)*sigma_e^16*gamma_s^7 + (1080*A^2*B^16*sigma_r^4 - 360*A*B^16*sigma_r^2 + 36*B^16)*sigma_e^16*gamma_s^6 + (5040*A^3*B^16*sigma_r^6 - 2160*A^2*B^16*sigma_r^4 + 360*A*B^16*sigma_r^2 - 24*B^16)*sigma_e^16*gamma_s^5 + (6*B^16 - 120*A*B^16*sigma_r^2 + 1080*A^2*B^16*sigma_r^4 - 5040*A^3*B^16*sigma_r^6 + 10080*A^4*B^16*sigma_r^8 + 10080*A^4*B^16*sigma_r^8*exp((gamma_s - 1)/(A*sigma_r^2)))*sigma_e^16*gamma_s^4 + 24*A*B^14*sigma_r^2*sigma_e^14*gamma_s^7 + (384*A^2*B^14*sigma_r^4 - 96*A*B^14*sigma_r^2)*sigma_e^14*gamma_s^6 + (2592*A^3*B^14*sigma_r^6 - 1152*A^2*B^14*sigma_r^4 + 144*A*B^14*sigma_r^2)*sigma_e^14*gamma_s^5 + (8064*A^4*B^14*sigma_r^8 - 5184*A^3*B^14*sigma_r^6 + 1152*A^2*B^14*sigma_r^4 - 96*A*B^14*sigma_r^2)*sigma_e^14*gamma_s^4 + (24*A*B^14*sigma_r^2 - 384*A^2*B^14*sigma_r^4 + 2592*A^3*B^14*sigma_r^6 - 8064*A^4*B^14*sigma_r^8 + 8064*A^5*B^14*sigma_r^10 + 8064*A^5*B^14*sigma_r^10*exp((gamma_s - 1)/(A*sigma_r^2)))*sigma_e^14*gamma_s^3 + 36*A^2*B^12*sigma_r^4*sigma_e^12*gamma_s^6 + (432*A^3*B^12*sigma_r^6 - 144*A^2*B^12*sigma_r^4)*sigma_e^12*gamma_s^5 + (2016*A^4*B^12*sigma_r^8 - 1296*A^3*B^12*sigma_r^6 + 216*A^2*B^12*sigma_r^4)*sigma_e^12*gamma_s^4 + (4032*A^5*B^12*sigma_r^10 - 4032*A^4*B^12*sigma_r^8 + 1296*A^3*B^12*sigma_r^6 - 144*A^2*B^12*sigma_r^4)*sigma_e^12*gamma_s^3 + (36*A^2*B^12*sigma_r^4 - 432*A^3*B^12*sigma_r^6 + 2016*A^4*B^12*sigma_r^8 - 4032*A^5*B^12*sigma_r^10 + 4032*A^6*B^12*sigma_r^12 + 4032*A^6*B^12*sigma_r^12*exp((gamma_s - 1)/(A*sigma_r^2)))*sigma_e^12*gamma_s^2 + 24*A^3*B^10*sigma_r^6*sigma_e^10*gamma_s^5 + (192*A^4*B^10*sigma_r^8 - 96*A^3*B^10*sigma_r^6)*sigma_e^10*gamma_s^4 + (576*A^5*B^10*sigma_r^10 - 576*A^4*B^10*sigma_r^8 + 144*A^3*B^10*sigma_r^6)*sigma_e^10*gamma_s^3 + (1152*A^6*B^10*sigma_r^12 - 1152*A^5*B^10*sigma_r^10 + 576*A^4*B^10*sigma_r^8 - 96*A^3*B^10*sigma_r^6)*sigma_e^10*gamma_s^2 + (24*A^3*B^10*sigma_r^6 - 192*A^4*B^10*sigma_r^8 + 576*A^5*B^10*sigma_r^10 - 1152*A^6*B^10*sigma_r^12 + 1152*A^7*B^10*sigma_r^14 + 1152*A^7*B^10*sigma_r^14*exp((gamma_s - 1)/(A*sigma_r^2)))*sigma_e^10*gamma_s + 6*A^4*B^8*sigma_r^8*sigma_e^8*gamma_s^4 + (24*A^5*B^8*sigma_r^10 - 24*A^4*B^8*sigma_r^8)*sigma_e^8*gamma_s^3 + (72*A^6*B^8*sigma_r^12 - 72*A^5*B^8*sigma_r^10 + 36*A^4*B^8*sigma_r^8)*sigma_e^8*gamma_s^2 + (144*A^7*B^8*sigma_r^14 - 144*A^6*B^8*sigma_r^12 + 72*A^5*B^8*sigma_r^10 - 24*A^4*B^8*sigma_r^8)*sigma_e^8*gamma_s + (6*A^4*B^8*sigma_r^8 - 24*A^5*B^8*sigma_r^10 + 72*A^6*B^8*sigma_r^12 - 144*A^7*B^8*sigma_r^14 + 144*A^8*B^8*sigma_r^16 + 144*A^8*B^8*sigma_r^16*exp((gamma_s - 1)/(A*sigma_r^2)))*sigma_e^8)/(B^16*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^16*gamma_s^8 + 8*A*B^14*sigma_r^2*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^14*gamma_s^7 + 28*A^2*B^12*sigma_r^4*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^12*gamma_s^6 + 56*A^3*B^10*sigma_r^6*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^10*gamma_s^5 + 70*A^4*B^8*sigma_r^8*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^8*gamma_s^4 + 56*A^5*B^6*sigma_r^10*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^6*gamma_s^3 + 28*A^6*B^4*sigma_r^12*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^4*gamma_s^2 + 8*A^7*B^2*sigma_r^14*exp((gamma_s - 1)/(A*sigma_r^2))*sigma_e^2*gamma_s + A^8*sigma_r^16*exp((gamma_s - 1)/(A*sigma_r^2)))

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!