Solve the integral in the point where the function is singular

15 次查看（过去 30 天）
Dimani4
Dimani4 2023-8-14
评论： Dimani4 2023-8-17
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 个评论
Torsten
Torsten 2023-8-14
编辑：Torsten 2023-8-14
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
Dimani4 2023-8-14
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.

请先登录，再进行评论。

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

回答（2 个）

Torsten
Torsten 2023-8-14
编辑：Torsten 2023-8-14
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 个评论
显示 6更早的评论
Dimani4
Dimani4 2023-8-15
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
Torsten 2023-8-15
Ran in:
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.

请先登录，再进行评论。

David Goodmanson
David Goodmanson 2023-8-16
编辑：David Goodmanson 2023-8-16
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 个评论
显示 1更早的评论
David Goodmanson
David Goodmanson 2023-8-17
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
Dimani4 2023-8-17
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品


版本

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by