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

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

该函数是以下形式的误差函数：

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

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

我知道D 应该在 1e-11 左右 and 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 的起始猜测，否则拟合将不起作用。

似乎有两个问题：

1. 您的函数并未针对每个实际值进行定义。假设您只对正值感兴趣D and 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 https://mathworks.com/help/symbolic/matlabfunction.html将其转换为普通的 matlab 函数，或者只是将结果复制并粘贴到脚本中，这样您就不必在每次迭代中以符号方式计算雅可比。

   您应该设置该选项'SpecifyObjectiveGradient',true https://mathworks.com/help/optim/ug/lsqcurvefit.html#input_argument_d0e68725告诉 matlab 使用你的精确雅可比矩阵。

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

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