Использование lsqcurvefit от Matlab для расчета бесконечного релаксационного спектра

Question

Я пытаюсь повторить анализ, выполненный в публикации, и у меня возникли проблемы. Речь идет о попытке вычислить бесконечный спектр релаксации по реологическим данным.

Учитывая экспериментальный ряд данных (выше), авторы извлекают спектр релаксации (H (тау)), решая для разности суммы квадратов между экспериментальными (выше) и модельными данными (ниже).

Эта релаксация следует распределению, описываемому следующим соотношением.

Этот расчет был выполнен с помощью функции Matlabs lsqcurefit. Однако lsqcurvefit не принимает функцию в качестве разрешимого параметра. Я хотел бы знать, как это извлечение может быть выполнено с помощью этой функции Matlab.

РЕДАКТИРОВАТЬ: После комментариев Bentoys, вот разъяснение моего вопроса.

Я подробно опишу следующее только для G '' (Gdp) для экономии места.

У меня есть экспериментальные данные в следующем виде: вектор значений частоты (омега), xdata и вектор значений отклика Gdp (Gdp).

Я хотел бы рассчитать H (тау), и для этого мне понадобятся параметры, содержащиеся в этой функции. Это дает следующее выражение, которое мне нужно решить:

ВВП - это функция от омеги, а мои xdata - это вектор значений омеги, но я интегрирую по отношению к ln (тау). Кажется, это может быть возможной проблемой?

Кроме того, у меня нет четкого идеала ожидаемых начальных значений для 6 переменных, только для результирующего H (tau), поэтому я выбрал произвольные значения для начала. Я могу оптимизировать их относительные значения, если рассчитанные значения могут быть получены.

Из ваших предложений мой код Matlab выглядит следующим образом:

w = numdata(:,1); %w is omega (experimental xdata)
GdpExp = numdata(:,3); % response values (ydata)
x0 = [10,10,0.1,0.1,1,1]; % arbitrary intial values
H = @(x, xdata) x(1)*exp(-(xdata-log(x(3))).^2/x(5)^2/2)...
    + x(2)*exp(-(xdata-log(x(4))).^2/x(6)^2/2);

Gdp = @(A_1, A_2, tau_1, tau_2, sigma_1, sigma_2, w) ...
    integral(H([A_1, A_2, tau_1, tau_2, sigma_1, sigma_2], ...
    u).* w.*exp(u)./(1+w.^2.*exp(2*u)), -Inf, Inf);

lsqcurvefit(Gdp, x0, w, GdpExp);

Эта валюта дает следующую ошибку:

>> lsqcurvefit(Gdp, x0, w, GdpExp);
Not enough input arguments.

Error in

Inf_Spec_Test>@(A_1,A_2,tau_1,tau_2,sigma_1,sigma_2,w)integral
(H([A_1,A_2,tau_1,tau_2,sigma_1,sigma_2],u).*w.*exp(u)./(1+w.^2.*exp(2*u)),
-Inf,Inf)

Error in lsqcurvefit (line 202)
        initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot 
continue.

Правильно ли я считаю, что каждая функция является автономной, или имена A_1, A_2 и т. Д. Должны совпадать с x (1), x (2), или A_1 и т. Д. Просто обозначают для обозначения расчетные значения?

Экспериментальные данные и рассчитанный H (тау) должны напоминать следующие цифры.

Я также обнаружил следующее соотношение, которое можно использовать для преобразования вектора омега в вектор тау, что может помочь преодолеть несоответствие между текущими xdata и интегралом.

Answer 1

Предполагая, что у вас есть экспериментальные данные H(tau) в форме вектора Y, вы хотите соответствовать вашей модели H:

H = @(A_1, A_2, tau_1, tau_2, sigma_1, sigma_2, tau) ...
      A1*exp(-(log(tau)-log(tau_1)).^2/sigma_1^2/2) + A2*exp(-(log(tau)-log(tau_2)).^2/sigma_2^2/2);

lsqcurvefit ожидает функцию в виде f(x, xdata). Здесь xdata - это «данные времени» с именами tau и x всех коэффициентов вашей функции. Итак, функция перезаписи H дает:

H = @(x, xdata) x(1)*exp(-(log(xdata)-log(x(3))).^2/x(5)^2/2) + x(2)*exp(-(log(xdata)-log(x(4))).^2/x(6)^2/2);

Тогда я полагаю, что у вас есть первоначальное предположение для значений всех шести коэффициентов. Если нет, фитинг может не сходиться к соответствующему решению. Я называю это первоначальное предположение x0, которое является вектором длины 6. Чтобы получить соответствие, просто позвоните:

lsqcurvefit(H, x0, tau, Y)

Возможно, вы захотите добавить дополнительные параметры (например, границы параметров), см. Документацию lsqcurvefit для получения более подробной информации.

EDIT

Следуя подробностям в комментариях ниже, у вас есть вектор значений G+G' для нескольких значений omega (в векторе W).

Итак, функция, которую вы хотите установить, немного сложнее:

H = @(x, xdata) x(1)*exp(-(xdata-log(x(3))).^2/x(5)^2/2) + x(2)*exp(-(xdata-log(x(4))).^2/x(6)^2/2);

GGp = @(A_1, A_2, tau_1, tau_2, sigma_1, sigma_2, w) ...
      integral(H([A_1, A_2, tau_1, tau_2, sigma_1, sigma_2], u).* w.*exp(u).*(1+w.*exp(u))./(1+w.^2.*exp(2*u)), -Inf, Inf);

после небольшого изменения переменной u = ln(tau).

Тогда это точно то же самое, за исключением того, что теперь у вас есть другая подходящая функция:

lsqcurvefit(GGp, x0, W, Y)