Nonlinear Curve fitting with integrals

Question

0 个投票

I encountered a nonlinear fitting problem, and the fitting formula is shown in Equation (1), which includes two infinite integrals (in practice, the integration range can be set from 0.01E-6 to 200E-6).

In these formulas, except for x and y being vectors, all other variables are scalars, and Rmedian and sigma are the parameters to be fitted.

I found a related post and tried to write the code based on it, but it keeps reporting errors. The error message seems to indicate that the vector dimensions are inconsistent, preventing the operation from proceeding. However, these functions are all calculations for individual scalars.

Error using  / 
Matrix dimensions must agree.
Error in dsdmain>@(r)1/(2*r*sigma*sqrt(2*pi))*exp(-(log(2*r)-log(2*Rmean)^2)/(2*sigma^2)) (line 13)
    gauss = @(r) 1/(2*r*sigma*sqrt(2*pi))* exp( -(log(2*r)-log(2*Rmean)^2)/(2*sigma^2) );

My question is: Can I refer to the content of this post to solve my problem? If yes, what does this error message mean? If not, how should I resolve my problem? (Note: The range of Rmedian is 1E-6 to 5E-6)

After modifying the code according to Walter Roberson and Torsten's suggestion, the program no longer throws errors. But no matter how I set the initial values on my end, it always prompts:

    Initial point is a local minimum.
Optimization completed because the size of the gradient at the initial point 
is less than the default value of the optimality tolerance.
<stopping criteria details>
theta = initial point
   1.0e-05 *
    0.2000
    0.1000
Optimization completed: The final point is the initial point.
The first-order optimality measure, 0.000000e+00, is less than
options.OptimalityTolerance = 1.000000e-06.
Optimization Metric                                       Options
relative first-order optimality =   0.00e+00  OptimalityTolerance =   1e-06 (default)

I have checked all the formulas and the units of the variables, and I didn't find any problems.

-------------------------------------- Beblow is my code for issue reproduction ----------------------------------------------------------

function testmain()
clc
    function Kvec = model(param,xdata)
        % Vector of Kd for every xdata:
        Kvec = zeros(size(xdata));
        Rmean = param(1);
        Rstd = param(2);
        
        for i = 1:length(xdata)
            model = @(r) unified(xdata(i),r,Rmean,Rstd,delta,Delta,D,lambda);
            Kvec(i) = integral(model,0.01E-6,200E-6);
        end
    end
    function s = unified(g,R,Rmean,Rstd,delta,Delta,D,lambda)
        %unified Unified fitting model for DSD
        %   exponentional combined
        
        factor = 1./(2*R*Rstd*sqrt(2*pi)) *2 ; % int(P(r)) = 0.5,1/0.5=2
        p1 = -(log(2*R)-log(2*Rmean)).^2/(2*Rstd^2);
        c = -2*2.675E8^2*g.^2/D;
        
        tmp = 0;
        for il = 1:length(lambda)
            a2 = (lambda(il)./R).^2;
            an4 = (R/lambda(il)).^4;
            Psi = 2+exp(-a2*D*(Delta-delta))-2*exp(-a2*D*delta)-2*exp(-a2*D*Delta)+exp(-a2*D*(Delta+delta));
            tmp = tmp+an4./(a2.*R.*R-2).*(2*delta-Psi./(a2*D));
        end
        p2 = c*tmp;
        
        s = factor.*exp(p1+p2);
    end
Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
g = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
xdata = g;
ydata = [100, 91.16805426, 80.52955192, 67.97705378, 55.1009735,41.87307917, 30.39638776, 21.13515607, 13.7125649, 8.33083767, 5.146756077, 2.79768609];
ydata = ydata/ydata(1); % normalize
% Inital guess for parameters:
Rmean0 = 2E-6;
Rstd0 = 1E-6;
p0 = [Rmean0;Rstd0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
theta = lsqcurvefit(@model,p0,xdata,ydata,[0.01E-6;0.1E-6],[10E-6,2E-6])
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Follow Question

Answer 1

Torsten 2025-4-12

编辑：Torsten 2025-4-13

在 MATLAB Online 中打开

0 个投票

This code works for your supplied data:

Delta = 0.075;

delta = 0.002;

D = 0.098E-9;

gamma = 2.675E8;

Rmin = 0.01e-6;

Rmax = 200e-6;

npoints = 100000;

lambda = [2.0816 5.9404 9.2058 12.4044 15.5792 18.7426 21.8997 25.0528 28.2034 31.3521];

xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];

ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;

% Optimize

mu0 = (Rmin + Rmax)/2;

sigma0 = 0.2;

p0 = [mu0,sigma0];

format long

p = lsqcurvefit(@(p,xdata)fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints),p0,xdata,ydata,[Rmin,0],[Rmax,Inf])

Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.

p = 1×2

0.000002705618161 0.001475482837009

% Plot results

y = fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints);

hold on

plot(xdata,ydata)

plot(xdata,y)

hold off

grid on

%

function y = fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints)

mu = p(1);

sigma = p(2);

y = zeros(size(xdata));

n1 = ceil(npoints*(p(1)-Rmin)/(Rmax-Rmin));

n2 = npoints - n1;

r1 = linspace(Rmin,mu,n1);

r2 = linspace(mu,Rmax,n2);

r = [r1,r2(2:end)];

for i = 1:numel(xdata)

F = fun(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);

y(i) = trapz(r,F);

end

%

function I = fun(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)

alpha = lambda./r.';

Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...

2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...

exp(-alpha.^2*D*(Delta+delta));

lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2.*(r.').^2-2).*(2*delta-Psi./(alpha.^2*D)),2);

E = exp(lnE);

P = 1./(2*r.'*sigma*sqrt(2*pi)).*exp(-log(r.'/mu).^2/(2*sigma^2));

I = (2*E.*P).';

end

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2025-4-16

编辑：Torsten 2025-4-16

在 MATLAB Online 中打开

I modified the code you provided by replacing trapz with Gauss-Legendre integration, but the results turned out worse. Why is this happening?

Available codes try to use a small number of grid points together with a high-order method to get a good approximation of an integral. This assumes that the grid points are placed at the correct positions in the integration interval. But this so called "adaptive gridding" is difficult in case of your function: for small values of sigma, the interesting part of your function is concentrated in a very narrow region. Without further help for the integrator (e.g. by prescribing a very fine grid (as I did in trapz) or by defining Rmean as a Waypoint (as I tested with "integral" and which worked, too), it is very difficult for the integrator to find this region and will fail.

So summarizing: in your case, the integration method doesn't matter, but it's essential that enough grid points are placed at the relevant positions. This failed for MATLAB's "integral" without further support, and it will most probably fail with every other integrator as well without giving it a helping hand (surprisingly, "quadgk" worked even without the waypoints option).

Here is the "Waypoints" code for the use of "integral":

Delta = 0.075;

delta = 0.002;

D = 0.098E-9;

gamma = 2.675E8;

Rmin = 0.01e-6;

Rmax = 200e-6;

npoints = 100000;

lambda = [2.0816 5.9404 9.2058 12.4044 15.5792 18.7426 21.8997 25.0528 28.2034 31.3521];

xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];

ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;

% Optimize

mu0 = (Rmin + Rmax)/2;

sigma0 = 0.2;

p0 = [mu0,sigma0];

format long

p = lsqcurvefit(@(p,xdata)fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata,[Rmin,0],[Rmax,Inf])

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.1e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.8e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 7.6e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.3e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.6e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.9e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.5e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.3e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-10. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.1e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.8e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.

p = 1×2

0.000002705623324 0.002304341838077

% Plot results

y = fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda);

hold on

plot(xdata,ydata)

plot(xdata,y)

hold off

grid on

%

function y = fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)

mu = p(1);

sigma = p(2);

for i = 1:numel(xdata)

F = @(r)fun(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);

y(i) = integral(F,Rmin,Rmax,'Waypoints',mu);

end

%

function I = fun(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)

alpha = lambda./r.';

Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...

2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...

exp(-alpha.^2*D*(Delta+delta));

lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2.*(r.').^2-2).*(2*delta-Psi./(alpha.^2*D)),2);

E = exp(lnE);

P = 1./(2*r.'*sigma*sqrt(2*pi)).*exp(-log(r.'/mu).^2/(2*sigma^2));

I = (2*E.*P).';

end

Shannon 2025-4-17

The code I'm using adopted from P313 of https://porousmedia.rice.edu/resources/PhD%20Thesis.pdf , which implements 24-point Gauss-Legendre integration. But my data isn't producing correct results.

Shannon 2025-4-17

@Sam Chak How amazing, the Gauss–Kronrod quadrature actually works!

请先登录，再进行评论。

Answer 2

Walter Roberson 2025-4-3

在 MATLAB Online 中打开

1 个投票

integral() always passes in a vector of locations to calculate at, unless you use the 'ArrayValued', true option.

sig = integral(f,0.01E-6,200E-6)/integral(gauss,0.01E-6,200E-6);

refers to

f = @(r) gauss(r)*kernel(r,delta,Delta,D)*2;

with the input r being a vector, gauss( r) can be expected to be a vector and kernal(...) can be expected to be a vector, so you have set up for two vectors to be "*" together. As the two vectors are likely the same size, you can expect problems from the * operator which does algebraic matrix multiplication for which the "inner dimensions" must match. You need to use the .* operator.

Meanwhile,

    gauss = @(r) 1/(2*r*sigma*sqrt(2*pi))* exp( -(log(2*r)-log(2*Rmean)^2)/(2*sigma^2) );

with the input r being a vector, the 1/(2*r*sigma*sqrt(2*pi)) component is 1/(vector) which is going to fail because the / operator is the matrix-right-divide operator in which the number of columns on the two sides must match. You need to use the ./ operator.

Then the exp() is going to be a vector and you have a * operator again, and that's going to fail because the inner dimensions must be the same size... again you need to be using the .* operator instead of the * operator.

15 个评论
显示 13更早的评论隐藏 13更早的评论

Torsten 2025-4-7

编辑：Torsten 2025-4-7

It doesn't help if your function is different from 0 in only one point - the numerical value of the integral will remain 0. If you really expect "sigma" in the order of magnitude of 1e-6, then it is impossible to use a numerical integrator for your task. If you try to plot P(r), you will see that you almost get the dirac distribution that is Inf at r = Rmean and 0 else.

Go one step back and try plotting P(r)*E(r,x,delta,s,D) over r for some of your x-values and see up to which value of "sigma" you get a reasonable function that you think MATLAB could successfully integrate. Then see if this value is still acceptable for your fitting task to start with.

Just a sidenote: Are you sure about the 2 in the denominator of the prefactor for P(R) ? The integral evaluates to 0.5 whereas I think it should evaluate to 1 as the probability density function for a lognormal distribution.

On Mathematics Stack Exchange you write

αm is the positive root of the equation J3/2(αmR)=αmRJ5/2(αmR)

So you mean that for each value of R, you get a different vector αm with the first 10 positive roots of the above equation ?

Maybe you could give a link to a publication where the above model has been used successfully to fit the parameters sigma and Rmedian.

Torsten 2025-4-8

编辑：Torsten 2025-4-8

Yes, that's why I substituted am with R via am = lambda / R

Ah, you solve for the first 10 positive roots for lambda of the equation J_3/2(lambda) = lambda*J_5/2(lambda) .

Is y(i)/y(1) = P(D>=x(i)) where D is the droplet size in mu-m ?

Shannon 2025-4-9

R is the droplet size (radius) and the distribution P® is supposed to be a lognornal distribution. Eq(2)

请先登录，再进行评论。

Answer 3

Torsten 2025-4-8

编辑：Torsten 2025-4-9

在 MATLAB Online 中打开

0 个投票

It takes too long to run the code online - try if it produces reasonable results. Maybe you have to set lower and upper bounds for the parameters in "lsqcurvefit".

Note that "mu" and "sigma" are not the expected value and the standard deviation of the Lognormal Distribution. Further I think that "sigma" is dimensionless.

Delta = 0.075;

delta = 0.002;

D = 0.098E-9;

gamma = 2.675E8;

Rmin = 0.01e-6;

Rmax = 200e-6;

lambda = [2.0816 5.9404 9.2058 12.4044 15.5792 18.7426 21.8997 25.0528 28.2034 31.3521];

xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];

ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;

mu0 = (Rmin + Rmax)/2;

sigma0 = 0.2;

% Plot integrand for initial values for mu and sigma

N = 100;

r = linspace(Rmin,Rmax,N);

I = arrayfun(@(x)fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu0,sigma0),xdata,'UniformOutput',0);

plot(r,cell2mat(I))

grid on

% Optimize mu and sigma

p0 = [mu0,sigma0];

p = lsqcurvefit(@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata)

function y = fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)

mu = p(1);

sigma = p(2);

y = zeros(numel(xdata),1);

for i = 1:numel(xdata)

F = @(r) fun_int(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);

y(i) = integral(F,Rmin,Rmax);

end

function I = fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)

I = zeros(numel(r),1);

for i = 1:numel(r)

alpha = lambda/r(i);

Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...

2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...

exp(-alpha.^2*D*(Delta+delta));

lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2).*(2*delta-Psi./(alpha.^2*D)));

E = exp(lnE);

P = 1./(2*r(i)*sigma*sqrt(2*pi))*exp(-log(r(i)/mu).^2/(2*sigma^2));

I(i) = 2*E*P;

end

7 个评论
显示 5更早的评论隐藏 5更早的评论

Walter Roberson 2025-4-9

在 MATLAB Online 中打开

Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
gamma = 2.675E8;
Rmin = 0.01e-6;
Rmax = 200e-6;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;
mu0 = (Rmin + Rmax)/2;
sigma0 = 0.2;
% Plot integrand for initial values for mu and sigma
N = 100;
r = linspace(Rmin,Rmax,N);
I = arrayfun(@(x)fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu0,sigma0),xdata,'UniformOutput',0);
whos r I
size(I{1})
plot(r,cell2mat(I(:)))
grid on
% Optimize mu and sigma
p0 = [mu0,sigma0];
p = lsqcurvefit(@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata)
function y = fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)
  mu = p(1);
  sigma = p(2);
  y = zeros('like',xdata);
  for i = 1:numel(xdata)
    F = @(r) fun_int(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);
    y(i) = integral(F,Rmin,Rmax);
  end
end
function I = fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)
  I = zeros('like', r);
  for i = 1:numel(r)
    alpha = lambda/r(i);
    Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...
          2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...
          exp(-alpha.^2*D*(Delta+delta));
    lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2).*(2*delta-Psi./(alpha.^2*D)));
    E = exp(lnE);
    P = 1./(2*r(i)*sigma*sqrt(2*pi))*exp(-log(r(i)/mu).^2/(2*sigma^2));
    I(i) = 2*E*P;
  end
end

However, running this generates a long list of messages similar to

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 9.3e+45. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

> In integralCalc>iterateScalarValued (line 393)

In

integralCalc>vadapt (line 148)

In

integralCalc (line 77)

In

integral (line 87)

In

q2175948>fun (line 32)

In

q2175948>@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda) (line 24)

In

lsqcurvefit/objective (line 314)

In finDiffEvalAndChkErr

In finitedifferences

In computeFinDiffGradAndJac

In

sfdnls (line 51)

In

snls (line 345)

In

lsqncommon (line 178)

In

lsqcurvefit (line 299)

In

q2175948 (line 24)

Notice how large the error bound is. There is no hope of fixing error bounds that large with small tweaks.

Shannon 2025-4-10

Origin can perform the fitting, but MATLAB fails? ‌

From one of the authors:

We did fit the dataset you provided to test in the way presented in the paper. Assuming a monomodal log-normal distribution (volume based) and using 10 roots (and not including free component or other effects), we get RMedian = 2.706 μm +/- 0.031 μm and log normal width factor of 0.0024 +/- 2.8231

The integration was included in the Origin Fitting function and defined using a LabTalk Script but see the help file in Origin help for integration notes. If you go to Origin help – All Books – Tutorials – Data Analysis (i.e., 4 Data Analysis) and look in the section on nonlinear fitting there are some sheets there that should help you create your function in Origin. Look at Fitting with Integral using LabTalk Function and Fitting with Two Integrals using LabTalk Function. There are some other files or ways in there, that may also help.
The help indicates that Origin also allows for +inf/-inf for limits, however, as mentioned we limited our integral as discussed after consideration of the probed range. Depending on the limits you use and the integral functions you may also need to split the range as indicated in the help file.
The origin help files we have indicate that using integrals needs version 8.6SR0 (or 8.0SR6 depending on the help file page) and up.
So, we defined the two integrals and then called on them in the function and instead of using +inf we use an upper limit as mentioned. The integrals themselves are done using Origin capability (numerically) (see also the help file for the math functions in LabTalk for integral and other files there). The function summation (with roots) was written in expanded form that way the roots were entered (fixed parameters) up to the number used but there is likely other ways as well.

Shannon 2025-4-11

I send my data in this post to the authors and they give me the fitting result back. The uncertainty in the width-factor is strange, but at least it can complete the fitting, and the fitted y matches the original y-vector very well.

Shannon 2025-4-11

xdata ydata [%] Fitting [%]

0.300616 100 100

0.53884 91.16805426 91.83196659

0.771392 80.52955192 80.6936292

1.009616 67.97705378 67.4689156

1.24784 55.1009735 53.86768785

1.480392 41.87307917 41.40418259

1.718616 30.39638776 30.28783311

1.95684 21.13515607 21.24404364

2.189392 13.7125649 14.45360593

2.427616 8.33083767 9.377799625

2.66584 5.146756077 5.865511151

2.898392 2.79768609 3.587601962

请先登录，再进行评论。

Nonlinear Curve fitting with integrals

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（2 个）

15 个评论
显示 13更早的评论隐藏 13更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论

类别

标签

Community Treasure Hunt

Nonlinear Curve fitting with integrals

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论 隐藏 3更早的评论

更多回答（2 个）

15 个评论 显示 13更早的评论 隐藏 13更早的评论

7 个评论 显示 5更早的评论 隐藏 5更早的评论

类别

标签

另请参阅

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论