required assistance in fitting

HI

I was trying to fit my data (in attchment) with fiist column V and second column J(V) with following equation (which honestly found complex to me)

mb = 9.11E-31; e = 1.6E-19; t = 3E-9; ℏ = 1.05E-34;

with Φ1 and Φ2 as the fitting parameters

Your support in this regard is apreiable!

Best

Somnath

[Mathieu NOE]

do we have an initial guess or range for Φ1 and Φ2 ?

[Somnath Kale]

@Mathieu NOE values of Φ1 and Φ2 are within the range of 0.1 to 5 in eV units

Answer 1

[Mathieu NOE]

so this is the poor"s man optimization (w/o any toolbox)

a brute force approach , let's create a 2 parameter grid and evaluate the function vs your data

the complex equation is simplified by taking some large blocks and define intermediate variables , instead of writing a unreadable long stuff.

we get a 2 dimensionnal error map (plotted after log conversion for better rendering) , and we pick the point of minimal error

we are fortunate that the error surface is rather smooth (banana)

for my own fun I did 2 iterations, the second is just refined grid around the first obtained minima point

% load J and V data
data = readmatrix('data.xlsx'); % first column V and second column J(V) 
V = data(:,1);
J = data(:,2);


mb = 9.11E-31;
e = 1.6E-19;
t = 3E-9;
h = 1.05E-34;

A  = 4*e*mb/(9*pi^2*h^3);

const_alpha = 4*t*sqrt(2*mb)/(3*h);

%Φ1 and Φ2 are within the range of 0.1 to 5 in eV units

%% First iteration
increment1 = 0.017;
increment2 = 0.019;
phi1 = (-min(0.5*V)+eps:increment1:2);
phi2 = (max(0.5*V):increment2:4);

for ci = 1:numel(phi1)
    for cj = 1:numel(phi2)
        
        P1 = phi1(ci)*e;
        P2 = phi2(cj)*e;
        Jm = model(const_alpha,e,A,P1,P2,V);
        err2(ci,cj) = mean((J - Jm).^2);
        
    end
end

[val,ind] = min(err2,[],'all','linear');
error_after_first_iteration = val

error_after_first_iteration = 3.9231e-07

range_factor = 100;
err2(err2>range_factor*val) = range_factor*val;

figure,imagesc(log10(err2));colorbar('vert');
colormap('jet')

% find optimum point
[r,c] = ind2sub(size(err2),ind);

P1 = phi1(r)*e;
P2 = phi2(c)*e;
Jm1 = model(const_alpha,e,A,P1,P2,V);

%% second iteration ?? 
clear err2
phi1 = linspace(phi1(r)*0.8,phi1(r)*1.2,100);
phi2 = linspace(phi2(c)*0.8,phi2(c)*1.2,100);

for ci = 1:numel(phi1)
    for cj = 1:numel(phi2)
        
        P1 = phi1(ci)*e;
        P2 = phi2(cj)*e;
        Jm = model(const_alpha,e,A,P1,P2,V);
        err2(ci,cj) = mean((J - Jm).^2);
        
    end
end

[val,ind] = min(err2,[],'all','linear');
error_after_second_iteration = val

error_after_second_iteration = 3.8205e-07

range_factor = 100;
err2(err2>range_factor*val) = range_factor*val;

figure,imagesc(log10(err2));colorbar('vert');
colormap('jet')

% find optimum point
[r,c] = ind2sub(size(err2),ind);

P1 = phi1(r)*e;
P2 = phi2(c)*e;
Jm2 = model(const_alpha,e,A,P1,P2,V);

figure,
plot(V,J,V,Jm1,V,Jm2)
xlabel('V');
ylabel('J(V)');
legend('raw','1st iteration','2nd iteration');

% optimal phi1 and phi2 (in eV)
phi1(r)

ans = 0.8005

phi2(c)

ans = 2.1355

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function 
function Jm = model(const_alpha,e,A,P1,P2,V)

        E = sqrt(P2 - 0.5*e*V);
        F = sqrt(P1 + 0.5*e*V);

        G = (P2 - 0.5*e*V).^1.5;
        H = (P1 + 0.5*e*V).^1.5;
        
        alpha = const_alpha./(P1 - P2 +e*V);

        N = exp(alpha.*(G - H));

        D = (alpha.^2).*(E - F).^2;

        B = 3/4*alpha.*e.*V;

        Jm = -A*N./D.*sinh(B.*(E-F));
end

[Somnath Kale]

Hi @Mathieu NOE

Thank you for responding! I'm still working on grasping the process. I'll reach out again once I've understood it better. Appreciate your help!

Regards, Somnath

[Mathieu NOE]

here, for my own fun, I used fminsearch (no toolbox required !) in the second iteration

still we need the first iteration, based on grid search principle , to find good IC for the fminsearch final optimization

again, similar results as above , the second iteration loop brings no big improvement, but it shows how to use fminsearch in this context

have fun !

error_after_first_iteration = 3.9231e-07

P1 = 0.7650

P2 = 2.1840

error_after_second_iteration = 3.8202e-07

P1 = 0.7999

P2 = 2.1360

% load J and V data
data = readmatrix('data.xlsx'); % first column V and second column J(V) 
V = data(:,1);
J = data(:,2);

mb = 9.11E-31;
e = 1.6E-19;
t = 3E-9;
h = 1.05E-34;

A  = 4*e*mb/(9*pi^2*h^3);

const_alpha = 4*t*sqrt(2*mb)/(3*h);

%Φ1 and Φ2 are within the range of 0.1 to 5 in eV units

%% First iteration
increment1 = 0.017;
increment2 = 0.019;
phi1 = (-min(0.5*V)+eps:increment1:2);
phi2 = (max(0.5*V):increment2:4);

for ci = 1:numel(phi1)
    for cj = 1:numel(phi2)
        
        P1 = phi1(ci)*e;
        P2 = phi2(cj)*e;
        Jm = model(const_alpha,e,A,P1,P2,V);
        err2(ci,cj) = mean((J - Jm).^2);
        
    end
end

[val,ind] = min(err2,[],'all','linear');
error_after_first_iteration = val

range_factor = 100;
err2(err2>range_factor*val) = range_factor*val;

figure,imagesc(log10(err2));colorbar('vert');
colormap('jet')

% find optimum point
[r,c] = ind2sub(size(err2),ind);

% optimal phi1 and phi2 (in eV)
P1 = phi1(r)
P2 = phi2(c)

Jm1 = model(const_alpha,e,A,P1*e,P2*e,V);

%% second iteration ?? 
% optimisation with fminsearch
% global const_alpha e A V J
IC = [phi1(r) phi2(c)]; 
[x,FVAL] = fminsearch(@objectivefcn1,IC);

% optimal phi1 and phi2 (in eV)
P1 = x(1)
P2 = x(2)
Jm2 = model(const_alpha,e,A,P1*e,P2*e,V);
error_after_second_iteration = FVAL

figure,
plot(V,J,V,Jm1,V,Jm2)
xlabel('V');
ylabel('J(V)');
legend('raw','1st iteration','2nd iteration');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function 
function Jm = model(const_alpha,e,A,P1,P2,V)

        E = sqrt(P2 - 0.5*e*V);
        F = sqrt(P1 + 0.5*e*V);

        G = (P2 - 0.5*e*V).^1.5;
        H = (P1 + 0.5*e*V).^1.5;
        
        alpha = const_alpha./(P1 - P2 +e*V);

        N = exp(alpha.*(G - H));

        D = (alpha.^2).*(E - F).^2;

        B = 3/4*alpha.*e.*V;

        Jm = -A*N./D.*sinh(B.*(E-F));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = objectivefcn1(x)
global const_alpha e A V J

    % P1 = phi1(r)*e;
    % P2 = phi2(c)*e;
    P1 = x(1)*e;
    P2 = x(2)*e;

    Jm = model(const_alpha,e,A,P1,P2,V);
    f = mean((J - Jm).^2);

end

[Mathieu NOE]

hello

problem solved ?