Numerically solve of integral equation with known Integrand

I'm trying to solve the attached integral equation (Eq. 5), numerically. The value of particle acceleration and fluid acceleration are known up to the time step n where t = n and the integral is from time 0 to t. I'm currently using the basic definition of integral to solve it (namely, Eq. 4).

The integral can be evaluated numerically except at a region where τ is approaching t. There, the integrand becomes infinite and needs to be treated separately. Therefore, I split the integral into two terms: Z1 and Z2 where Z1 is the series term expressed in Eq. 5 and Z2 is the last term in Eq. 5. Here is the code that I wrote:

if n > 2
    Z1= 0;
    for m = 2:n-1
         Z1 = Z1+ ((d2z_p(n-m) - d2z_f(n-m))/sqrt(m)) * dt^(1/2)
    end
    Z2 = (d2z_p(n-1) - d2z_f(n-1)) * dt^(1/2);
    Z_1 = Z1 + Z2;
end

Z_1 is supposed to be the solution of the integral, but I know that it doesn't give the right answer either because Eq. 4 is not the right one to use, or my code is wrong. Anyone can help me on how to solve the given integral?

Answer 1

[Torsten]

If d2z_p(i) = dv_p / dtau at tau = i * t / n = i * dt and d2z_f(i) = dv_f / dtau at tau = i * t / n = i * dt in your code from above, it is correct.

Note that you can also use

Z= 0;
for m = 1:n-1
   Z = Z+ ((d2z_p(n-m) - d2z_f(n-m))/sqrt(m)) * dt^(1/2)
end

to approximate the integral without splitting in Z1 and Z2.

[AHMAD TALAEI]

Thank you, Torsten. At least I know the code is correct. Do you know if there is another way to solve that integral?

[Torsten]

t = ...;      % vector of t-values where dvpdt and dvfdt are available
dvpdt = ...;  % dvpdt vector with the same length as t vector 
dvfdt = ...;  % dvfdt vector with the same length as t vector
t1 = ...;     % upper limit of the integral
t0 = 0.0;     % lower limit of the integral
eps = 1.0e-6;
fun=@(tau)1./sqrt(t1-tau).*(interp1(t,dvpdt,tau)-interp1(t,dvfdt,tau));
sol = integral(fun,t0,t1-eps)

[AHMAD TALAEI]

Thanka a lot!

Both give the same answer, although the "integral" function slows down the code. I just noticed that I have multiplied the integral by a wrong coefficient and that is the reason it didn't give the correct solution.

Ahmad