Integral Recurrence Formula Implementation

10 次查看（过去 30 天）
Djordje
Djordje 2025-4-18
编辑： Torsten 2025-4-18
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
  3 个评论
显示 1更早的评论
Djordje
Djordje 2025-4-18
I provided more context by editing the original post.
Torsten
Torsten 2025-4-18
编辑：Torsten 2025-4-18
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" ?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2025-4-18
编辑：Torsten 2025-4-18
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by