finding fit parameters between two functions

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)

[Torsten]

In my opinion, it's a very good accordance between measurements and simulation - especially the times when the peaks occur. Was the blue curve achieved without parameter fitting ? I wouldn't go this step of further adaption.

[Eli]

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]

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

Answer 1

[Torsten]

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.

[Eli]

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?

[Torsten]

Your function must accept the array of values for tau for which you have measurement data "times_tau" and the parameter vector "parameters" as input and should give the modelled data "modelled_data_at_times_tau" as output:

function [modelled_data_at_times_tau] = untitled5(parameters,times_tau)

and the call to lsqcurvefit should be

parameters = lsqcurvefit(@untitled5, parameter0, times_tau, measured_data_at_times_tau);

This is explained in detail if you look at the examples provided under

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

[Eli]

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?

[Torsten]

Are you sure a variable in the order of 1e-16 influences "interpolated_IM_FINAL" ?

You could test it by computing the numerical derivative of "interpolated_IM_FINAL" with repect to x(5) for the initial parameter vector.

For this, call "untitled5" with x(5) and x(5)+1e-4 and compute the numerical gradient

(interpolated_IM_FINAL([1600,1600,0.5,14000,0.1+1e-4,4]) - interpolated_IM_FINAL([1600,1600,0.5,14000,0.1,4]))/1e-4

[Eli]

why in the order of 1e-16?

In the function, I'm changing

IM_Final=E_tau(IM_Final,time_shift);

Where the time_shift is my variable x0(5).

you mean because of the multiplication of 1e-15? when I'm using the internal function E_tau, it's shifting the function in the time domain and as a result the 'interpolated_IM_FINAL' and 'average_data_set' will be closer or farther apart

1e-4 might not change, but 0.2 and 0.1 is noticably different.

Note that I also tried with 0.002 as x(5) just in case a step size is 1e-4, and than multiplied x(5)*1e-15*100 just in case, and still the same, x(5) doens't change.

[Torsten]

What do you get as dfj/dx(5), the derivative (i.e. sensitivity) of "interpolated_IM_FINAL" with respect to x(5) for the initial parameter vector ? I told you above how to approximately compute it.

Depending on whether the vector,

v = (interpolated_IM_FINAL([1600,1600,0.5,14000,0.1+1e-4,4]) - interpolated_IM_FINAL([1600,1600,0.5,14000,0.1,4]))/1e-4    

consists of small or large elements, the change in x(5) will be small or large. My guess is that v = 0.

The 1e-4 is related to the options setting "FiniteDifferenceStepSize". You must choose a larger value for the FiniteDifferenceStepSize for x(5) than the default (default is sqrt(eps), I guess).

sqrt(eps)      

ans = 1.4901e-08

[Eli]

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]

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

Glad to hear that it seems to work now.