finding fit parameters between two functions

Question

0 个投票

Hello

I have a script, where i have a simulation that creates a function.

On the other hand I have a measurement. They don't fit and I want to use fminsearch to try to find the variables that can make them fit.

The idea is to find the variables N_set , A_dir, A_seq and maybe some offset in time (x axis)and Const (y axis).

Would appriciate some help.

The simulation and data are presented. this is example of one specific value set, where we can see the don't have a very good overlap.

Would appriciate your help.

Thanks.

load('matlab1.mat')

untitled5()

Unrecognized function or variable 'time_average_data'.

Error in untitled5 (line 169)
hold on; plot(time_average_data,average_data_set-min(average_data_set),'LineWidth', 1.5)

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

Eli 2023-12-23

编辑：Eli 2023-12-23

So I just played with the parameters by hand. that's why i need parameter optimization.

Without this, it's totally off so more precise fit will bring me closer to my real parameters. (will have many similar experiments so optimization should help a lot)

Image Analyst 2023-12-23

Is the function that gives a fit just a single mathematical formula? If so, what is it?

请先登录，再进行评论。

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

Follow Question

Answer 1

Torsten 2023-12-23

编辑：Torsten 2023-12-23

2 个投票

Make a function that - given tau and the parameters you want to fit - returns I_m(tau).

Then look at the examples for "lsqcurvefit" to see how you can use this function as the "fun" in the call to the optimizer.

https://uk.mathworks.com/help/optim/ug/lsqcurvefit.html

Use the parameters you fitted by hand as vector of initial values for the parameter vector.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Eli 2023-12-23

编辑：Eli 2023-12-23

在 MATLAB Online 中打开

Trying this but I'm doing something wrong.

this is my function:

function [interpolated_IM_FINAL,time_array,IM_Final] = untitled5(x0)
A_dir = x0(1);
A_seq = x0(2);
N_set = x0(3);
Const = x0(4);

I'm interested in finding those optimal values:

A_dir,A_seq,N_set,Const

Trying something like this:

clear
clc
load('matlab1.mat');
x0 = [1600,1600,0.5,14000]; %initial guess
x = lsqcurvefit(untitled5,x0,time_average_data,average_data_set)

or

x = lsqcurvefit(@untitled5, x0, time_average_data, average_data_set)

None works.

Q: What am I doing wrong here?

Google bard "told" me i need to define

anonymous_function = @(x0, time_average_data, average_data_set) untitled5(x0);

and than use:

x = lsqcurvefit(anonymous_function, x0, time_average_data, average_data_set);

Q: why this is what i need to do?

after running what bard suggested, i get:

"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.

<stopping criteria details>"

Q: what does that mean?

Eli 2023-12-23

编辑：Torsten 2023-12-23

在 MATLAB Online 中打开

seems to work now- thanks.

Just one question left:

I'm attaching now the new function and running the following code:

load('matlab1.mat');
x0 =  [1600,1600,0.5,14000,0.1,4];
x = lsqcurvefit(@untitled5, x0, time_average_data, average_data_set);
% Clear variables and command window
%function [interpolated_IM_FINAL,IM_Final,time_array] = untitled5(x0,time_average_data)
function interpolated_IM_FINAL = untitled5(x0,time_average_data)
% A_dir = 1600;
% A_seq = 1600;
% Const = 200;
% N_set = 0.6;
% Constants and Parameters
%A_dir,A_seq,N_set,Const
A_dir = x0(1);
A_seq = x0(2);
N_set = x0(3);
Const = x0(4);
time_shift = x0(5)*1e-15;
pulse_time = x0(6);
% time_average_data = 2*1e15*(80e-3*25*10e-9/3e8)*[-49:50];
c = 299792458; % Speed of light [m/s]
twidth = 20e-15; % Time window [s]
n = 2^13; % Number of points
lambda0 = 756*1e-9; % Center wavelength [m]
W0 = 2*pi*c/lambda0; % Center frequency [rad/s]
FWHM = pulse_time*1e-15; % Initial pulse duration [s]
Epsilon_0 = 8.8541878128*1e-12;
% Time domain
T = linspace(-twidth/2, twidth/2, n); % Time grid
dt = T(2) - T(1);
Dnu = 1/dt;
delta_T = twidth/n; % Time interval [s]
% Gaussian pulses in time domain
sigma = FWHM/((2*log(2))^0.5);
gaussian = exp(-T.^2/sigma^2);
% Time difference between small pulses (time period of the wave)
pulse_period = lambda0 / c;
t_shift = pulse_period / 2;
FMHM_new = 250e-18;
sigma_new = FMHM_new/((2*log(2))^0.5);
% Calculate the height of the original Gaussian at the centers of the new Gaussians
height_gaussian_center = exp(0); % Gaussian is maximum at its center
% Normalize the pulses based on the original Gaussian envelope
N1 = height_gaussian_center / max(exp(-T.^2/sigma^2));
% Calculate the value of the original Gaussian at the position of gaussian2
value_at_gaussian2 = exp((t_shift)^2/sigma^2);
N2 = height_gaussian_center / value_at_gaussian2;
N2= N_set;
% Calculate the value of the original Gaussian at the position of gaussian3
%value_at_gaussian3 = exp((t_shift)^2/sigma^2);
%N3 = height_gaussian_center / value_at_gaussian3;
N3 = N2;
gaussian1 = exp(-T.^2/sigma_new^2) * N1;
gaussian2 = exp(-(T - t_shift).^2 / sigma_new^2) * N2;
gaussian3 = exp(-(T + t_shift).^2 / sigma_new^2) * N3;
xuv1_lamda = 62e-9; % nm   %pump
W0_xuv1 = 2*pi*c/xuv1_lamda; % Center frequency [rad/s]
xuv2_lamda = 27.5e-9; % nm  %probe
W0_xuv2 = 2*pi*c/xuv2_lamda; % Center frequency [rad/s]
f_in_t = gaussian1 + gaussian2 + gaussian3;
% Combine three Gaussians
fin_t_1 = gaussian1 .* exp(1i*W0_xuv1*T);
fin_t_2 = gaussian2 .* exp(1i*W0_xuv1*(T - t_shift));
fin_t_3 = gaussian3 .* exp(1i*W0_xuv1*(T + t_shift));
fin_t_xuv1 = fin_t_1 + fin_t_2 + fin_t_3;
fin_t_1 = gaussian1 .* exp(1i*W0_xuv2*T);
fin_t_2 = gaussian2 .* exp(1i*W0_xuv2*(T - t_shift));
fin_t_3 = gaussian3 .* exp(1i*W0_xuv2*(T + t_shift));
fin_t_xuv2 = fin_t_1 + fin_t_2 + fin_t_3;
% Function to calculate E(t-tau) using circshift
E_tau = @(E, tau) circshift(E, round(tau/dt));
% Initialize array for IM(tau) and use T directly as tau values
IM_tau_array_seq = zeros(size(T));
IM_tau_array_direct = zeros(size(T));
% Perform numerical integration for IM(tau)
%figure;
% Original Gaussian at 756nm
% subplot(4, 1, 1);
% plot(1e15*T, gaussian, 1e15*T,-gaussian,'LineWidth', 1.5,'Color',"blue");
% hold on
% plot(1e15*T, real(gaussian.*exp(i*W0*T)), 'LineWidth', 1.5);
% xlabel('Time (fs)');
% ylabel('Amplitude');
% title('NIR field at 756nm');
% grid minor
% 
% % XUV Pulses
% subplot(4, 1, 2);
% plot(1e15*T, abs(fin_t_xuv1),1e15*T, -abs(fin_t_xuv1), 'LineWidth', 1.5,'color',"blue");
% hold on 
% plot(1e15*T, real(fin_t_xuv1), 'LineWidth', 1.5);
% xlabel('Time (s)');
% ylabel('Amplitude');
% title(['XUV Pulse ' num2str(1240/(xuv1_lamda*1e9)) 'eV']);
% grid minor
% 
% subplot(4, 1, 3);
% plot(1e15*T, abs(fin_t_xuv2),1e15*T, -abs(fin_t_xuv2), 'LineWidth', 1.5,'color',"blue");
% hold on 
% plot(1e15*T, real(fin_t_xuv2), 'LineWidth', 1.5);
% xlabel('Time (s)');
% ylabel('Amplitude');
% title(['XUV Pulse ' num2str(round(1240/(xuv2_lamda*1e9))) 'eV']);
% grid minor
% IM(tau)
%subplot(4, 1, 4);
for i = round(length(T)/8):round(length(T) - length(T)/8)
    weight = zeros(1,length(T));
    tau = T(i);
% Initialize IM(tau) for the current tau
    
    % Calculate E(t-tau) using circshift
    E_t_tau1 = E_tau(fin_t_xuv1, -tau);
    E_t_tau2 = E_tau(fin_t_xuv2, -tau);
    % Calculate IM(tau)
    IM_t_tau2 = abs(E_t_tau2).^2; 
    IM_t_tau1 = abs(E_t_tau1).^2; 
    weight(1)=trapz(T(1:1),IM_t_tau1(1:1)+(1:1),2)*dt;
    % for j=2:length(T)
    % % Update IM(tau)
    % weight(j)=trapz(T(1:j),IM_t_tau1(1:j)+IM_t_tau2(1:j))*dt;
    % end
    cumulative_sum = cumtrapz(T, IM_t_tau1 + IM_t_tau2) * dt;
    weight(2:end) = cumulative_sum(2:end);
    IM_tau_array_seq(i) = trapz(T, abs(fin_t_xuv2).^2.*weight) * dt;
    %IM_t_tau = abs((fin_t_xuv1 + E_t_tau2).^2).^2;
    IM_t_tau = abs(abs(fin_t_xuv1).^2 + abs(E_t_tau2)).^2;
    IM_tau_array_direct(i) = trapz(T, IM_t_tau) * dt;
end
 IM_tau_array_seq(1:round(length(T)/8)-1) = IM_tau_array_seq(round(length(T)/8));
 IM_tau_array_seq(round(length(T) - length(T)/8)+1:end) = IM_tau_array_seq(round(length(T) - length(T)/8));
 IM_tau_array_direct(1:round(length(T)/8)-1) = IM_tau_array_direct(round(length(T)/8));
 IM_tau_array_direct(round(length(T) - length(T)/8)+1:end) = IM_tau_array_direct(round(length(T) - length(T)/8));
IM_tau_array_direct = (IM_tau_array_direct-min(IM_tau_array_direct))/max((IM_tau_array_direct-min(IM_tau_array_direct)));
%plot(1e15*T, IM_tau_array_direct, 'LineWidth', 1.5);
 IM_tau_array_seq = (IM_tau_array_seq - min(IM_tau_array_seq))/max((IM_tau_array_seq - min(IM_tau_array_seq)));
% figure();
% plot(1e15*T, A_dir.*IM_tau_array_direct+A_seq.*IM_tau_array_seq+Const, 'LineWidth', 1.5);
% xlabel('\tau');
% ylabel('I_M(\tau)');
% title('S (\tau) for all \tau');
% hold on; plot(time_average_data,average_data_set-min(average_data_set),'LineWidth', 1.5)
% grid minor
IM_Final=A_dir.*IM_tau_array_direct+A_seq.*IM_tau_array_seq+Const;
time_array = T.*1e15;
IM_Final=E_tau(IM_Final,time_shift);
interpolated_IM_FINAL = interp1(time_array, IM_Final, time_average_data);
end

But for some reason x(5) doesn't seem to be optimized at all (doesn't chnage it's value). don't understand why. could you please take a look?

Eli 2023-12-24

编辑：Eli 2023-12-24

在 MATLAB Online 中打开

Perfect!

That solved it.

The name by the way is " FinDiffRelStep" instead of "FiniteDifferenceStepSize"

I used than (didn't find how to only set one variable X(5)) :

options = optimset('FinDiffRelStep', [sqrt(eps) sqrt(eps) sqrt(eps) sqrt(eps) 1e-2 sqrt(eps)]);
x = lsqcurvefit(@untitled5, x0, time_average_data, average_data_set,lb,ub,options);

And, it seems to start changing. thanks.

Torsten 2023-12-24

If you use "optimoptions" instead of "optimset", the name is "FiniteDifferenceStepSize".

Glad to hear that it seems to work now.

请先登录，再进行评论。

finding fit parameters between two functions

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

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

类别

标签

Community Treasure Hunt

finding fit parameters between two functions

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

采纳的回答

8 个评论 显示 6更早的评论 隐藏 6更早的评论

更多回答（0 个）

类别

标签

另请参阅

Community Treasure Hunt

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

8 个评论
显示 6更早的评论隐藏 6更早的评论