Assistance with integration limits and implementing it in Matlab

Greetings all,

   Back at it again with a little more complicated issue that I tried, but stuck in neutral here.

   I'm going to introduce the functions first, and then how I coded them, and then ask my question because, while it makes sense to me equation wise, coding it is another thing (per usual)!

   So I have the following integral function (sorry I can't do symbols here):

    Integral with limits Psi_n-1 (read "Psi sub n-1) to infinity Phi(Psi)dPsi - Integral with limits Psi_n to infinity Phi(Psi)dPsi

   where Phi(Psi) = 4*Psi/(4(Psi)^2)^3/2

   and Psi = z/r. 

    Now for how I coded this:

      r=3.66;
       z=[1.3336 1.3370 1.3438 1.3594 1.3700 1.3577 1.3507 1.3559 1.2900 1.3525];

       Psi=z/r;

        n=10;

       Phi=@(Psi) (4*Psi_/(4*(Psi).^2+1).^3/2;
       q_Phi_a=integral(Phi,n-1,inf, 'ArrayValued', true);
       q_Phi_b=integral(Phi,n,inf, 'ArrayValued', true);

   Now here's the thing. q_Phi_a and q_Phi_b appear to be the same thing, so when you subtract them, there going to be zero, which is not what I thought. The only thing I can think of is how I am doing the limits.

    In the equation, the bottom limit is either Psi sub n-1 or Psi sub n. I don't think I am doing this right is Matlab. Any advice on this? Do I need another function handle in the limits of the integral function?

     Let me know if I need to attach anything from MathType for a clearer picture.

    Thanks!
    -J

Answer 1

[Torsten]

As written, your integral is simply

1/(4*Psi(n-1)^2+1)^0.5 - 1/(4*Psi(n)^2+1)^0.5

(an antiderivative of f(x)=4x/(4x^2+1)^1.5 is F(x)=-1/(4x^2+1)^0.5)

Best wishes

Torsten.

[Jesse]

Torsten,

Thanks for the response. But what about the limits of integration?

[Jesse]

Ah, ok I see what's going on, but what about the fact that Psi changes with respect to z? I mean Psi is a vector, so if I set this equation equal to q, q should also be a vector.

[Torsten]

The above value of the integral is

Integral with limits Psi_n-1 (read "Psi sub n-1) to infinity Phi(Psi)dPsi - Integral with limits Psi_n to infinity Phi(Psi)dPsi.

Make a loop over n to get the vector you are searching for:

for n=2:numel(Psi)
value_integral(n-1)=1/(4*Psi(n-1)^2+1)^0.5 - 1/(4*Psi(n)^2+1)^0.5;
end

Best wishes

Torsten.