# Solve the integral in the point where the function is singular

12 visualizzazioni (ultimi 30 giorni)
Dimani4 il 14 Ago 2023
Commentato: Dimani4 il 17 Ago 2023
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.  ##### 2 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Dimani4 il 14 Ago 2023
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.

Accedi per commentare.

### Risposte (2)

Torsten il 14 Ago 2023
Modificato: Torsten il 14 Ago 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 CommentiMostra 7 commenti meno recentiNascondi 7 commenti meno recenti
Torsten il 15 Ago 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))
value_integral =
0.069558476982290537581125337200736
So you see: at least for this case, your approach is better than mine.

Accedi per commentare.

David Goodmanson il 16 Ago 2023
Modificato: David Goodmanson il 16 Ago 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 CommentiMostra 2 commenti meno recentiNascondi 2 commenti meno recenti
Dimani4 il 17 Ago 2023
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.

Accedi per commentare.

### Categorie

Scopri di più su Numerical Integration and Differentiation in Help Center e File Exchange

R2016b

### Community Treasure Hunt

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

Start Hunting!