Integral Recurrence Formula Implementation

I have integral

where q is positive even number and C1>0 is a constant. This integral has recurrcive formula solution in the form:

with the base case:

I created this matlab function to implement this recurrence formula, but it seems not to work when compared against numerical integration. Here is my code and any help is very appriciated. One hint is that when q becomes large (in my case it can get to 50 to 60), the numbers become very large and I am not sure if my reccurence formula loses precision due to Matlab rounding large numbers. If necessary, I can provide bigger context for this code and how this function fits into bigger picture. Thank you in advance.

Code context

C1 is 100 x 100 array of positive numbers (values around between 100 and 1000), Y is an array (100 x 100 in size) with values between -5000 and 5000. S = 2500. Further context is that I have a cirle of radius S and my lower limits of the integral are points on the circomference of the circle in the 2nd and 3rd quadrant, while xP are all points in the circle.

function Iq = compute_Iq(C1, q, xP, Y, S)
% Function to compute Iq using recurrence formula

% Precompute common terms
CS = (xP.^2 + Y.^2) <= S^2;
xE = -sqrt(S^2 - Y.^2);
xE(~CS) = 0;
xP(~CS) = 0;
XC1E = (xE.^2 + C1.^2).^(3/2);
XC1P = (xP.^2 + C1.^2).^(3/2);

% Recurrence iterations
if q > 0
    c = 1; T = cell(1,q/2); C = T;
    for i = q:-2:2
        % Evaluate the first term of each integration by parts at the limits
        T1xP = xP.^(i-1) .* XC1P;
        T1xE = xE.^(i-1) .* XC1E;
        T{c} = (1/(i+2)) * (T1xP - T1xE);

        % Coefficient of the second term of the integration by parts
        C{c} = (i-1)*C1.^2/(i+2);
        c = c+1;
    end
end

% Base cases for q = 0. q = 1 is not possible becasue q is always even.
% int sqrt(x^2 + C1^2) dx = (x/2) * sqrt(x^2 + C1^2) + (C1^2/2) * log(x + sqrt(x^2 + C1^2))
I0P = (xP *  sqrt(xP.^2 + C1.^2) + C1.^2 * log(xP + sqrt(xP.^2 + C1.^2)))/2;
I0E = (xE *  sqrt(xE.^2 + C1.^2) + C1.^2 * log(xE + sqrt(xE.^2 + C1.^2)))/2;
Iq = I0P - I0E;

% Compuite integral Iq by accumulating all recurrence terms
if q > 0
    N = length(T);
    for i = N:-1:1
        Iq = T{i} - C{i}.*Iq;
    end
end
end

[Torsten]

Could you explain the inputs Y and S of your function ?

[Djordje]

I provided more context by editing the original post.

[Torsten]

C1 is 100 x 100 array of positive numbers (values around between 100 and 1000), Y is an array (100 x 100 in size) with values between -5000 and 5000. S = 2500. Further context is that I have a cirle of radius S and my lower limits of the integral are points on the circomference of the circle in the 2nd and 3rd quadrant, while xP are all points in the circle.

I still don't understand what you try to compute for Iq from these inputs. What are upper and lower limits of the (100x100 ?) different integral values you try to compute ? Why do you say that the lower limits of the integral are "points" ?

Answer 1

[Torsten]

This seems to work:

xl = 0;
xu = 4;
q = 12;
C1 = 2;
vl = (xl*sqrt(xl^2+C1^2)+C1^2*log(xl+sqrt(xl^2+C1^2)))/2;
vu = (xu*sqrt(xu^2+C1^2)+C1^2*log(xu+sqrt(xu^2+C1^2)))/2;
qq = 2;
for i = 1:q/2
    vl = xl^(qq-1)*(xl^2+C1^2)^(3/2)/(qq+2)-(qq-1)*C1^2/(qq+2)*vl;
    vu = xu^(qq-1)*(xu^2+C1^2)^(3/2)/(qq+2)-(qq-1)*C1^2/(qq+2)*vu;
    qq = qq + 2;
end
format long
vu-vl

ans = 
     2.178774669469912e+07

syms x
f = x^q*sqrt(x^2+C1^2);
double(int(f,x,xl,xu))

ans = 
     2.178774669469912e+07