Solve the integral in the point where the function is singular

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.

[Torsten]

The question is: Does the multiplication with "y_appr" result in a function "func_approach" that is integrable between x_approach(1) and x_approach(end) ?

[Dimani4]

Hi Torsten,

One function (y_appr) I get from measurements so I cannot fit it to polynom and one function (func1) which known. So the solution can be numerical integral with trapz function but if you have singularity in your function you cannot do trapz because you get Infinite. So my solution was to divide the integral into two parts. One part is to solve it symblolically in the point of singularity and one part is to solve it numerically in the rest of the integral.

Here I attach actually what should I solve.

Maybe this attachment clears the things up.

Thank you.

Answer 1

[Torsten]

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

[Torsten]

Example:

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: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

value_integral = 

0.069591339330869360294116043473469

func = sqrt(0.025)/(8*sqrt(pi*(x-x_approach(1))))*x^2;
vpaintegral(func,x,0,3)

ans = 

0.0695294

[Dimani4]

Thank you Torsten.

Can you please explain in more details what you did? I understood you do linear fit every time for two adjacent points and then you solve the integral. But why you didnt make ordinary fit for 2 points like:

Piecewise linear fit needs if you have points scattering and when just linear fit is not good enough then you divide your x for segments and you do for every segment linear fit. But here you are dealing only with 2 adjacent points.

Thank you very much.

[Torsten]

Instead of using "trapz", I did the same as you did for the first interval [x_approach(1),x_approach(2)] also for the subsequent points [x_approach(i),x_approach(i+1)] and added the results for the integrals.

"expr" is the equation of the line connecting (x_approach(i),y_approach(i)) and (x_approach(i+1),y_approach(i+1)), and I multiplied this line equation with your func1_symb and integrated it over [x_approach(i),x_approach(i+1)].

I don't understand what you mean by

I understood you do linear fit every time for two adjacent points and then you solve the integral. But why you didnt make ordinary fit for 2 points like:

What's the difference between the two options ?

[Dimani4]

Torsten,

I did with

expr=((y_approach(i+1)-y_approach(i))/(x_approach(i+1)-x_approach(i)))*x+(x_approach(i+1)*y_approach(i)-x_approach(i)*y_approach(i+1))/(x_approach(i+1)-x_approach(i));

As the regular linear approximation with two points and got exactly the same answer.

[Torsten]

Fine. I prefer the Lagrange polynomial representation because it generalizes for polynomial interpolation in more than two points:

https://en.wikipedia.org/wiki/Lagrange_polynomial

And it's easier to remember :-)

[Dimani4]

Torston,

With your approach I got a number:0.19552314913231515409487715206893

and with my previous approach to fit only the problematic point and then to do trapz I got: 0.19571766254784561667179332289379.

[Dimani4]

Torston,

You already answered to this question: one is Lagrange polynom (your approach) and the second one is the regular linear fit (mine).

"I understood you do linear fit every time for two adjacent points and then you solve the integral. But why you didnt make ordinary fit for 2 points like:

What's the difference between the two options ?"

Thank you very much.

[Torsten]

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))

value_integral = 

0.069558476982290537581125337200736

So you see: at least for this case, your approach is better than mine.

Answer 2

[David Goodmanson]

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.

[Dimani4]

Thank you for that approach. But I cannot see here how you bypass this problematic, singular point. You have "y" as Omega function of t in the mine integral which is measured, you cannot fit it. The approach of Torsten and me is to linear fit the function in the singular point and the point after and then multiply this linear function with 1/sqr(t-z). Then this integral can be solved analytically and then sum it with the rest which is solved with trapz function. My question is; is this is the right approach? Maybe there is some mathematically known approach for such a situation.

Thank you.

[David Goodmanson]

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.

[Dimani4]

Hi David,

func1=sqrt(0.025)./(8*sqrt(pi*(x_approach-x_approach(1))));

func2=(y1-y1(1)).*func1;

wholeintegral=trapz(x_approach(1:end),func2(1:end))

wholeintegral =

NaN

This is what you meant?

Thank you.