Matlab 曲线拟合不适用于小值（1e-12），我该怎么办？

Question

我安装了曲线拟合工具箱，我正在尝试将扩散数据拟合到特定函数。

函数是一个错误函数，形式如下：

y = 3500 - 2500 * erf( ( x-x0 ) / ( 2 * sqrt( D * t )) )

我希望应用程序为我提供 D 和 x0 的合理值，而 t 是预定义的常数。拟合所基于的数据点包括 x 和 y 的值。

我知道 D 应该在 1e-11 左右，而 x0 应该在 0.0014 左右，但该函数不会自行找到这些解决方案。尝试使用曲线拟合应用程序的默认参数时，Matlab 输出错误“输入必须真实且完整。”。 当我将 x0 的起始猜测设置为 0.0014 或 0.0015 时，它会找到正确的解决方案。但只有这两个值。为了找到 D 的正确解，我需要在方程中设置一个前置因子，就像在这个例子中一样 (1e-12)：

y = erf( ( x-x0 ) / ( 2 * sqrt( 1e-12 * D * t )) )

通过这种方式，matlab 可以找到正确的解决方案，但仅适用于 1e-10 到 1e-13 之间的前置因子。

这是非常有问题的，因为 D 的正确解决方案将在 1e-3 和 1e-15 之间变化，具体取决于我要使用的数据集。 x0 的值也会有所不同。所以以这种方式我无法实现一个通用的解决方案。

您对如何处理这个问题有什么建议吗？我不敢相信matlab不能解决这个问题，一定有办法。是不是因为数值太小了？

这是我正在使用的示例数据集：

y = [6000 6000 6000 6000 6000 6000 6000 6000 6000 5750 5500 5250 5000 4500 4000 3250 2750 2250 1750 1500 1400 1250 1250 1150];
x = [0:0.0001:0.0023];

当在方程中使用以下固定参数时，结果线非常适合数据点。但 matlab 不会找到它们。

D = 7.1e-11;
t = 900;
x0 = 0.0015;

（请记住，这些参数基于比我在此处提供的更大、更准确的数据集）

任何帮助都会很棒！非常感谢你。

这是一个工作示例，其中包含所有预定义参数以使适合工作。代码由曲线拟合工具（cftool）生成：

%% Fit: 'untitled fit 1'.

   % Input data
    C = [6000 6000 6000 6000 6000 6000 6000 6000 6000 5750 5500 5250 5000 4500 4000 3250 2750 2250 1750 1500 1400 1250 1250 1150]';
    x = [0:0.0001:0.0023]';

    [xData, yData] = prepareCurveData( x, C );

    % Set up fittype and options.
    ft = fittype( '3500-2500*erf((x-x0)/(2*sqrt(1e-10*D*900)))', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.MaxIter = 4000;
    opts.StartPoint = [0.5 0.0014];
    opts.Upper = [1 Inf];

    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts )

    % Plot fit with data.
    figure( 'Name', 'untitled fit 1' );
    h = plot( fitresult, xData, yData );
    legend( h, 'C vs. x', 'untitled fit 1', 'Location', 'NorthEast' );
    % Label axes
    xlabel x
    ylabel C
    grid on

（仅为方便起见绘制）

请注意，我在 sqrt() 项中包含了前因子 1e-10，并且我使用 0.0014 作为 x0 的起始猜测，否则拟合将不起作用。

Answer 1

似乎有2个问题：

1. 您的函数并未针对每个实数值进行定义。假设您只对 D 和 x0 的正值感兴趣，您可以通过将下限指定为 [0 0] 来简单地将搜索范围绑定到正数

2. 由于值非常小，Matlab 在数值上评估函数的导数时存在问题。因此，最好的解决方案是象征性地计算函数的雅可比。这可以使用符号工具箱来完成：

   syms x x0 D t real;
   
   y = 3500 - 2500 * erf( ( x-x0 ) / ( 2 * sqrt( D * t )) )
   J = jacobian(y, [D, x0]);
   

   您可以使用matlabFunction 将其转换为普通的matlab 函数，或者只是将结果复制并粘贴到您的脚本中，这样您就不必在每次迭代中符号地计算雅可比。

   您应该设置选项 'SpecifyObjectiveGradient',true 来告诉 matlab 使用您的确切 jacobian。

---

实现这两种解决方案会产生以下代码：

ydata = [6000 6000 6000 6000 6000 6000 6000 6000 6000 5750 5500 5250 5000 4500 4000 3250 2750 2250 1750 1500 1400 1250 1250 1150];
xdata = 0:0.0001:0.0023;
t = 900;

D = 0.1;
x0 = 0.1;

options = optimoptions('lsqcurvefit', 'SpecifyObjectiveGradient',true);
X = lsqcurvefit(@(x, xdata) y(x(1), x(2), t, xdata),[D, x0], xdata, ydata, [0 0], [], options);

D = X(1); % 6.4833e-11
x0 = X(2); % 0.0015

figure
hold on
plot(xdata, ydata);
plot(xdata, y(D, x0, t, xdata));

function [F, J] = y(D, x0, t, x)
  F = 3500 - 2500 * erf( ( x-x0 ) / ( 2 * sqrt( D * t )) );
  J =  [(1250.*t.*exp(-(x - x0).^2/(4.*D.*t)).*(x - x0))/(pi^(1/2).*conj((D.*t).^(3/2))); ...
            (2500.*exp(-(x - x0).^2/(4.*D.*t)))/(pi.^(1/2).*conj((D.*t).^(1/2)))]';
end