Error in cumtrapz. Input arguments to function include colon operator. To input the colon character, use ':' instead.

Hi

I would like to use the cumulative trapezoidal method to obtain the integral of function f with a domain [0.2, 1] with regards to x. However, it is showing errors in cumtrapz.

Any help is much appreciated! Thank you

This is my code:

syms Swi kmin kmax D muyg muyw Krgwi sat x 
Swi = 0.2;
kmin = 0.1;
kmax = 1;
D = 1;
muyg = 0.1;
muyw = 1;
Krgwi = 0.7;
sat = linspace(Swi,1,100);
f =@(x)(0.5.*(kmin - kmax) + kmax - 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)).*(0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg) .* 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax -0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2.*kmax - (2.*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(3/2))./((Swi - x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)));
cumtrapz(f,sat)

UPDATED: Thanks for all your help.

I'm given the function f and the goal is to take the integral of f in the domain of [0.2, 1], denoted as ksi

my goal is to plot the result of ksi in terms of x, with x ranges from 0.2 to 1. Initially I thought of the cumtrapz as the function seems difficult to perform analytically. Please correct me if I'm wrong.

[Stephen23]

CUMTRAPZ performs numerical integration: its documentation specifies that its inputs must be numeric arrays.

What do you expect CUMTRAPZ to do with your symbolic inputs?

[Dyuman Joshi]

Why define the variables as symbolic variables when you are just over-writing them?

As @Stephen23 has clarified, cumtrapz() performs numerical integration, where as you have provided a function handle.

Edit - Can you specify what exactly are you trying to do? Can you mention the mathematical formulation of the problem you are trying to solve?

[Sam Chak]

Hi @Ivan

Is the function plotted correctly over the range for ? The surface for the bivariate function is also plotted. Can you verify that?

syms Swi kmin kmax D muyg muyw Krgwi sat x 
Swi   = 0.2;
kmin  = 0.1;
kmax  = 1;
D     = 1;
muyg  = 0.1;
muyw  = 1;
Krgwi = 0.7;
% sat   = linspace(Swi, 1, 100);
sat   = 1;
f     = @(x) ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*( (1 - x)./(1 - Swi) ).^2 + kmax*(1 - x)./(1 - Swi)).*( 0.5*(kmin - kmax)*( (1 - x)./(1 - Swi) ).^2 + kmax.*(1 - x)./(1-Swi) )./(((muyw*Krgwi*D./muyg) .* 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax -0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2.*kmax - (2.*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(3/2))./((Swi - x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi))./( ( (muyw*Krgwi*D/muyg)*0.5*(kmin - kmax) .* ( (1 - x)./(1 - Swi) ).^2 + kmax*(1 - x)./(1 - Swi) ) + ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi) ) )

f = function_handle with value:
    @(x)(0.5*(kmin-kmax)+kmax-0.5*(kmin-kmax)*((1-x)./(1-Swi)).^2+kmax*(1-x)./(1-Swi)).*(0.5*(kmin-kmax)*((1-x)./(1-Swi)).^2+kmax.*(1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*0.5.*(kmin-kmax).*((1-x)./(1-Swi)).^2+kmax.*(1-x)./(1-Swi))+(0.5.*(kmin-kmax)+kmax-0.5.*(kmin-kmax).*((1-x)./(1-Swi)).^2+kmax.*(1-x)./(1-Swi))).*(kmax-kmin)./((kmax+kmin).*(Swi-1)*((2.*kmax-(2.*(kmax-kmin).*(sat-1))./(Swi-1))./(kmax+kmin)).^(3/2))./((Swi-x)./(1-Swi))+(0.5.*(kmin-kmax)+kmax-0.5.*(kmin-kmax).*((1-x)./(1-Swi)).^2+kmax*(1-x)./(1-Swi))./(((muyw*Krgwi*D/muyg)*0.5*(kmin-kmax).*((1-x)./(1-Swi)).^2+kmax*(1-x)./(1-Swi))+(0.5*(kmin-kmax)+kmax-0.5*(kmin-kmax)*((1-x)./(1-Swi)).^2+kmax*(1-x)./(1-Swi)))

fplot(f, [0.2 1]), grid on
xlabel('x'), ylabel('f(x)')

[x, sat] = meshgrid(0.2:0.01:0.4);
f     = ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*( (1 - x)./(1 - Swi) ).^2 + kmax*(1 - x)./(1 - Swi)).*( 0.5*(kmin - kmax)*( (1 - x)./(1 - Swi) ).^2 + kmax.*(1 - x)./(1-Swi) )./(((muyw*Krgwi*D./muyg) .* 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax -0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2.*kmax - (2.*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(3/2))./((Swi - x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi))./( ( (muyw*Krgwi*D/muyg)*0.5*(kmin - kmax) .* ( (1 - x)./(1 - Swi) ).^2 + kmax*(1 - x)./(1 - Swi) ) + ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi) ) );
surf(x, sat, f, 'FaceAlpha', 0.5, 'EdgeColor', 'none')
xlabel('x'), ylabel('sat'), zlabel('f(x, sat)')

[Ivan]

Thanks @Dyuman @Stephen23. I had updated my goal in the post. I would like to perform the integral with respect to x and to plot the result of the integral versus x (from 0.2 to 1). Can I do it with cumtrapz or should I perform another method?

[Walter Roberson]

Your f(x) involves the term (sat-1) but sat is a vector so f(x) returns a vector for scalar x

[Sam Chak]

Hi @Ivan

First things first, could you mathematically verify if your function appears as shown below? You certainly don't want to integrate the wrong function. Since the variable "sat" is defined over the range [0.2, 1.0], I assume that "sat" is, in fact, the variable "x". Additionally, "sat" only appears once in the function, which might be an unintended typo. Please correct me if I'm mistaken.

syms Swi kmin kmax D muyg muyw Krgwi sat f(x)
Swi   = 0.2;
kmin  = 0.1;
kmax  = 1;
D     = 1;
muyg  = 0.1;
muyw  = 1;
Krgwi = 0.7;
% sat = linspace(Swi, 1, 100);
sat   = x;  % assumed that "sat" is, in fact, the variable "x"
f(x)  = (0.5.*(kmin - kmax) + kmax - 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)).*(0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg) .* 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax -0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2.*kmax - (2.*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(3/2))./((Swi - x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)))

f(x) = 

[Ivan]

@Sam Chak Thank you for pointing that out!! Indeed I made a type and sat is indeed variable x as I made a mistake in the previous step!

[Sam Chak]

Since you are working in the symbolic realm, try using int(). Let us know the result.

https://www.mathworks.com/help/symbolic/sym.int.html

[Walter Roberson]

Q = @(v) sym(v);
syms Swi kmin kmax D muyg muyw Krgwi sat f(x)
Swi   = Q(0.2);
kmin  = Q(0.1);
kmax  = Q(1);
D     = Q(1);
muyg  = Q(0.1);
muyw  = Q(1);
Krgwi = Q(0.7);
% sat = linspace(Swi, 1, 100);
sat   = x;  % assumed that "sat" is, in fact, the variable "x"
f  = (Q(0.5).*(kmin - kmax) + kmax - Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)).*(Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg) .* Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax - Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((Q(2).*kmax - (Q(2).*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(Q(3)/2))./((Swi - x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax-Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax-Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)))

f = 

T = linspace(Swi+eps, 1, 150);
F = double(subs(f, x, T));
plot(T, F)

limit(f, x, Swi)

ans = 

NaN

vpaintegral(f, x, Swi, 1, 'MaxFunctionCalls', 1e9)

Error using sym/vpaintegral
Failed precision goal. Try using 'MaxFunctionCalls'.

At the lower end of the range, there is a 0/0 in an area that is otherwise headed to -infinity. Numeric integration is not able to get a solution with a reasonable number of calls.

Answer 1

[Walter Roberson]

Change

cumtrapz(f,sat)

to

cumtrapz(sat, f(sat))

[Dyuman Joshi]

I tried using cumtrapz(f(sat)) while writing my comment to see what the result would be. I didn't include in my comment, but it gives the same output as cumptrapz(sat,f(sat))

Swi = 0.2;
kmin = 0.1;
kmax = 1;
D = 1;
muyg = 0.1;
muyw = 1;
Krgwi = 0.7;
sat = linspace(Swi,1,100);
f =@(x)(0.5.*(kmin - kmax) + kmax - 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)).*(0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg) .* 0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax -0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2.*kmax - (2.*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(3/2))./((Swi - x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (0.5.*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)));

out1=cumtrapz(sat,f(sat));
out2=cumtrapz(f(sat));

disp(out1')

     0
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf

disp(out2')

     0
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf
   Inf

[Walter Roberson]

It turns out that the expression for f(x) involves a term (sat-1) which is the entire sat array. That should probably be a reference to either x or Swi instead of sat.

[Ivan]

Thank you Walter, the term sat is indeed x. Sorry as I didn't double-check it sooner

[Walter Roberson]

By the way, if you simplify(expand(f)) then you get out a notably shorter expression . However it still has the same integration characteristics. Maple says that if you start from Swi that both versions are undefined. Probably because of the singularity at 19/27 - (2*sqrt(994))/135

Answer 2

[Sam Chak]

Hi @Ivan

Perhaps you can learn from this example to compute the cumulative integral of univariate function, .

x = linspace(-4, 4, 8001);
f = 1./cosh(x).^2;  % contains function values for f(x) = 1/cosh²(x) in the domain [-4, 4].
plot(x, f), grid on
xlabel('x'), ylabel('f(x)')
title({'$f(x) = \frac{1}{\cosh^{2}(x)}$'}, 'interpreter', 'latex', 'fontsize', 16)

g = cumtrapz(x, f);     % g(x) = tanh(x) + 1
plot(x, g), grid on
xlabel('x'), ylabel('g(x)')
title({'Cumulative distribution function of $f(x)$'}, 'interpreter', 'latex', 'fontsize', 16)

[Sam Chak]

Hi @Ivan

As you probably already know, the function has two singularities over [0.2, 1.0]. MATLAB's integral() can evaluate and return some value. However, I am unsure whether the integral is convergent or divergent.

https://www.mathworks.com/help/matlab/ref/integral.html

Swi   = 0.2;
kmin  = 0.1;
kmax  = 1;
D     = 1;
muyg  = 0.1;
muyw  = 1;
Krgwi = 0.7;
f     = @(x) ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*((1 - x)/(1 - Swi)).^2 + kmax*(1 - x)/(1 - Swi)).*(0.5*(kmin - kmax)*((1 - x)/(1 - Swi)).^2 + kmax*(1 - x)/(1 - Swi))./(( (muyw*Krgwi*D/muyg)*0.5*(kmin - kmax)*((1 - x)/(1 - Swi)).^2 + kmax*(1 - x)/(1 - Swi)) + (0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*((1 - x)/(1 - Swi)).^2 + kmax*(1 - x)/(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((2*kmax - (2*(kmax - kmin)*(x - 1))/(Swi - 1))/(kmax + kmin)).^(3/2))./((Swi - x)/(1 - Swi)) + (0.5*(kmin - kmax) + kmax-0.5.*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi))./( ( (muyw*Krgwi*D/muyg)*0.5*(kmin - kmax) .* ( (1 - x)./(1 - Swi) ).^2 + kmax*(1 - x)./(1 - Swi) ) + ( 0.5*(kmin - kmax) + kmax - 0.5*(kmin - kmax)*((1 - x)./(1 - Swi)).^2 + kmax*(1 - x)./(1 - Swi) ) );
q     = integral(f, 0.2, 1.0)

Warning: Minimum step size reached near x = 0.2. There may be a singularity, or the tolerances may be too tight for this problem.

q = -1.9936e+03

x     = linspace(0.2, 1.0, 801);
plot(x, f(x), 'linewidth', 1.5, 'Color', "#0437F2"), grid on
xlabel('x'), ylabel('f(x)')

[Sam Chak]

Hi @Ivan

The following approach employs numerical integration using vpintegral() as demonstrated by @Walter Roberson in the comment. What distinguishes it is that this variable precision-based integration evaluates from 0.2+ε to 1.0, yielding a finite value. Nevertheless, I have reservations about the accuracy of the result, given the presence of a singularity at .

Q     = @(v) sym(v);
syms Swi kmin kmax D muyg muyw Krgwi sat f(x)
Swi   = Q(0.2);
kmin  = Q(0.1);
kmax  = Q(1);
D     = Q(1);
muyg  = Q(0.1);
muyw  = Q(1);
Krgwi = Q(0.7);
% sat = linspace(Swi, 1, 100);
sat   = x;  % assumed that "sat" is, in fact, the variable "x"
f     = (Q(0.5).*(kmin - kmax) + kmax - Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)).*(Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg) .* Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax - Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))) .* (kmax - kmin)./((kmax + kmin).*(Swi - 1)*((Q(2).*kmax - (Q(2).*(kmax - kmin).*(sat - 1))./(Swi - 1))./(kmax + kmin)).^(Q(3)/2))./((Swi - x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax-Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi))./(((muyw*Krgwi*D./muyg).*Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)) + (Q(0.5).*(kmin - kmax) + kmax-Q(0.5).*(kmin - kmax) .* ((1 - x)./(1 - Swi)).^2 + kmax .* (1-x)./(1-Swi)))

f = 

% locating the singularities
den1  = 45/22*x - 5/22 == 0;
xs1   = vpasolve(den1)      % singularity point #1

xs1 = 

0.11111111111111111111111111111111

den2  = 5/4*x - 1/4 == 0;
xs2   = vpasolve(den2)      % singularity point #2

xs2 = 

0.2

den3  = 5/2*x + 27/10*(5/4*x - 5/4).^2 - 61/20 == 0;
xs3   = vpasolve(den3)      % singularity point #3 & #4

xs3 = 

% plot of f(x) showing the singularities at 4 locations
fplot(f, [0 1.2])
xlabel('x'), ylabel('f(x)')

% check if limits exist near the singularities xs2 and xs3
Lxs2  = limit(f, x, xs2+eps)

Lxs2 = 

Lxs3  = limit(f, x, xs3(1))

Lxs3 = 

% numerical integration from 0.2+eps to 1.0
vpaintegral(f, x, xs2+eps, 1, 'MaxFunctionCalls', 1e9)

ans = 

Looking at the previous approach, integral(f) returns .