Solve the integral in the point where the function is singular
18 views (last 30 days)
Show older comments
Hi people,
I need to solve some numerical integral where I have a singular point. Let me explain you what I mean:
I have some function which I measure (I cannot fit it to polynom), then I need to multiply this function to some another function which I know and then to do integral. So the function is :
func1=sqrt(0.025)./(8*sqrt(pi*(x_approach-x_approach(1))));
func_approach=func1.*y_appr;
So func1 I know and y_appr I have from measurements. I need to do inegral from minimum of x_approach to maximum of x_approach. My first idea was to do numerical integral everywhere except the point where the numerical integral becomes infinite (in the first point). The point where integral becomes infinite I thought to do regular symbolic integral.The function y_approach I thought to fit to polynom of first order and then multiply with the func1 and then to do inegral in the boundaries of from point1 (x_approach(1) where func1 is infinite) to point2 (x_approach(2)). Something like that:
p = polyfit(x_approach(1:2),y_appr(1:2),1);
syms x_appr
expr=p(1)*x_appr+p(2);
func1_symb=sqrt(0.025)/(8*sqrt(pi*(x_appr-x_approach(1))));
func_appr_symb=func1_symb*expr;
Fvpaint = vpaintegral(func_appr_symb,x_appr,[x_approach(1) x_approach(2)]);
Fvpaint_rest_ofintegral=trapz(x_approach(2:end),func_approach(2:end));
thewholeintegral=Fvpaint_rest_ofintegral+Fvpaint;
But as for me I have a doubt if my approach is correct or this approach is too much "rough". I guess I need to do this symbolic integral in the point of singularity in the boundaries of +-very small delta and then the rest just to solve numerical integral. The question how I do it and what delta should I take, how to estimate the delta?
Thank you very much.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456772/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456777/image.jpeg)
2 Comments
Answers (2)
Torsten
on 14 Aug 2023
Edited: Torsten
on 14 Aug 2023
I'd suggest
syms x
func1_symb=sqrt(0.025)/(8*sqrt(pi*(x-x_approach(1))));
value_integral = 0.0;
for i = 1:numel(x_approach)-1
expr = (x-x_approach(i))/(x_approach(i+1)-x_approach(i))*y_approach(i+1) + ...
(x-x_approach(i+1))/(x_approach(i)-x_approach(i+1))*y_approach(i); %piecewise linear
%expr=(y_approach(i)+y_approach(i+1))/2; % piecewise constant
func_appr_symb=func1_symb*expr;
Fvpaint = vpaintegral(func_appr_symb,x,[x_approach(i) x_approach(i+1)]);
value_integral = value_integral + Fvpaint;
end
value_integral
8 Comments
Torsten
on 15 Aug 2023
You could try "trapz" for the example above, compare the results and see which approach approximates the real value better:
x_approach = 0:0.1:3;
y_approach = x_approach.^2;
syms x
func1_symb=sqrt(0.025)/(8*sqrt(pi*(x-x_approach(1))));
value_integral = 0.0;
for i = 1:1
expr = (x-x_approach(i))/(x_approach(i+1)-x_approach(i))*y_approach(i+1) + ...
(x-x_approach(i+1))/(x_approach(i)-x_approach(i+1))*y_approach(i); %piecewise linear
%expr=(y_approach(i)+y_approach(i+1))/2; % piecewise constant
func_appr_symb=func1_symb*expr;
Fvpaint = vpaintegral(func_appr_symb,x,[x_approach(i) x_approach(i+1)]);
value_integral = value_integral + Fvpaint;
end
value_integral = value_integral + trapz(x_approach(2:end),sqrt(0.025)./(8*sqrt(pi*(x_approach(2:end)-x_approach(1)))).*y_approach(2:end))
So you see: at least for this case, your approach is better than mine.
David Goodmanson
on 16 Aug 2023
Edited: David Goodmanson
on 16 Aug 2023
Hi Dimani,
Leaving out the first '1 ' term that has no singularity and leaving off some constants, you are basically looking at
Int{z,b} y(t) / sqrt(t-z) dt
where I assume the lower limit is z every time for all choices of z. (I suppose the value of b may change as z changes but that does not affect the approach).
You have a set of t values and an experimental set of y values. I assume here that z is one of the values of t. If z falls in between two values of t the problem can still be done with interpolation, but there is an extra step involved.
If set of You can add and subtract y(z) / sqrt(t-z) to get
Int{z,b} (y(t)-y(z))/sqrt(t-z) dt + Int{z,b} y(z)/sqrt(t-z) dt
where the subtraction in the numerator of the first integral removes the singularity and the second integral is done analytically. Here us an example.
format compact
t = (0:.1:20);
zind = 51;
z = t(zind)
bind = 180;
b = t(bind)
y = 20 + 4*t; % example with an analytic solution for comparison
y1 = (y - y(zind))./sqrt(t-z);
% the value of y1 at t = z may come out inf or nan due to numerical
% precision issues or 0/0 form so make sure it's zero.
y1(zind) = 0;
I1 = trapz(t(zind:bind),y1(zind:bind))
I2 = 2*sqrt(b-z)*y(zind) % analytic
I3 = I1+I2
% analytic calculation of the whole thing, should agree with I3
I4 = 2*20*sqrt(b-z) + 4*(2/3)*(b-z)^(3/2) +4*z*2*sqrt(b-z)
(I4-I3)/I4
z = 5
b = 17.9000
I1 = 123.5272
I2 = 287.3326
I3 = 410.8597
I4 = 410.8856
ans = 6.2868e-05
The result is good to better than one part in 10^4 which seems reasonable for a calculation on only 201 points or less. The calculation could be much more accurate if you knew the derivative of y(t) at t=z but since your y is experimental, a truly precise value is not so easy to come by.
3 Comments
David Goodmanson
on 17 Aug 2023
I forgot to mention that I was assuming that z is one of the values of the t vector so I updated the answer I posted accordingly. So, with that assumption, then
y(t) / sqrt(z-t) has a singularity at t = z because the numerator y(z) is nonzero there. Now if y(z) / sqrt(t-z) is subtracted off, the resulting function (y(t) - y(z)) / sqrt(t-z) has numerator 0 at t = z, so there is no singularity and that function can be integrated numerically. Then y(z) sqrt(t-z) can be integrated analytically and the sum of those two integrals agrees with the exact answer by 6 parts in 10^5. The method definitely works. It's really the same method that you and Torsten are doing, but using just the consant term of the linear fit and applying it to the entire span of t (in a favorable situation for that) rather than a neighborhood of z.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!