How to fit my experimental data with using matlab code

I need to fit my experimental data ( α vs.ℏω) using following eq. (1) which includes an integral and several fitting parameters. I would be very appeciative if anyone can help me on this.

My experimental data is below. The first column is the photon energy ℏω, the second column is the absorption coefficient α.

3.50048 111558.1

3.49306 110412.3

3.48567 109296.9

3.4783 108206.6

3.47098 107140

3.46368 106095

3.45641 105068.7

3.44917 104061.7

3.44196 103072.1

3.43479 102099.8

3.42764 101143.8

3.42053 100205.1

3.41344 99281.64

3.40638 98374.32

3.39936 97482.5

3.39236 96606.17

3.38539 95745.09

3.37845 94898.81

3.37154 94066.66

3.36466 93249.77

3.3578 92447.08

3.35098 91658.52

3.34418 90883.3

3.33741 90122.05

3.33067 89373.87

3.32395 88639.14

3.31727 87917.02

3.31061 87207.67

3.30397 86510.62

3.29737 85825.79

3.29079 85152.95

3.28424 84492.14

3.27771 83842.59

3.27121 83204.22

3.26474 82577.2

3.25829 81961.28

3.25187 81355.93

3.24547 80761.22

3.2391 80176.73

3.23275 79603.08

3.22643 79039.7

3.22014 78486.55

3.21386 77945.1

3.20762 77414.7

3.2014 76897.3

3.1952 76394.22

3.18903 75908.27

3.18288 75443.4

3.17675 75004.61

3.17065 74598.71

3.16458 74235.05

3.15852 73923.8

3.15249 73678.54

3.14649 73513.53

3.1405 73443.66

3.13454 73483.41

3.1286 73644.72

3.12269 73934.94

3.1168 74354.01

3.11093 74892.9

3.10508 75529.38

3.09925 76229.77

3.09345 76945.23

3.08767 77614.51

3.08191 78166.34

3.07618 78523.23

3.07046 78607.37

3.06477 78345.64

3.05909 77675.59

3.05344 76552.35

3.04781 74951.98

3.0422 72870.59

3.03662 70331.23

3.03105 67377.13

3.0255 64069.29

3.01998 60486.95

3.01447 56713.59

3.00899 52842.91

3.00352 48961.84

2.99808 45156.05

2.99265 41497.55

2.98725 38048.77

2.98186 34854.71

2.9765 31949.08

2.97115 29346.73

2.96583 27053.16

2.96052 25059.89

2.95523 23350.06

2.94996 21901.12

2.94472 20684.62

2.93948 19671.4

2.93427 18831.74

2.92908 18137.02

2.92391 17561.13

2.91875 17080.8

2.91362 16676.22

2.9085 16330.32

2.9034 16029.45

2.89831 15762.78

2.89325 15521.56

2.88821 15299.38

2.88318 15091.21

2.87817 14893.54

2.87317 14703.65

2.8682 14519.72

2.86324 14340.42

2.8583 14164.89

2.85338 13992.52

2.84847 13822.86

2.84359 13655.74

2.83871 13490.86

2.83386 13328.17

2.82902 13167.64

2.8242 13009.14

2.8194 12852.66

2.81461 12698.16

2.80984 12545.68

2.80509 12395.1

2.80035 12246.43

[Torsten]

Use "lsqcurvefit" to fit your parameters and "integral" to approximate the integral ( given values of the parameters from "lsqcurvefit" ).

[Chaoping]

Dear Torsten,

Thank you for your kind reply. Indeed, I tried to use the "lsqcurvefit" to fit the parameters. I used following MATLAB

code (also see the attached file) for my problem. When I run this code, it appears "invalid use of operator on line

161". Would you please help me to take a look. Many thanks for your help.

Best regards,

Chaoping

function theta = Gil

data = [

3.50048 111558.0625

3.49306 110412.2578

3.48567 109296.8984

3.4783 108206.6484

3.47098 107140.0078

3.46368 106094.9922

3.45641 105068.6641

3.44917 104061.6953

3.44196 103072.125

3.43479 102099.7813

3.42764 101143.8281

3.42053 100205.0703

3.41344 99281.64063

3.40638 98374.32031

3.39936 97482.5

3.39236 96606.17188

3.38539 95745.09375

3.37845 94898.8125

3.37154 94066.66406

3.36466 93249.76563

3.3578 92447.07813

3.35098 91658.51563

3.34418 90883.30469

3.33741 90122.04688

3.33067 89373.86719

3.32395 88639.14063

3.31727 87917.01563

3.31061 87207.67188

3.30397 86510.61719

3.29737 85825.78906

3.29079 85152.94531

3.28424 84492.14063

3.27771 83842.59375

3.27121 83204.21875

3.26474 82577.19531

3.25829 81961.28125

3.25187 81355.92969

3.24547 80761.21875

3.2391 80176.72656

3.23275 79603.07813

3.22643 79039.69531

3.22014 78486.54688

3.21386 77945.10156

3.20762 77414.70313

3.2014 76897.30469

3.1952 76394.21875

3.18903 75908.26563

3.18288 75443.39844

3.17675 75004.60938

3.17065 74598.71094

3.16458 74235.05469

3.15852 73923.79688

3.15249 73678.53906

3.14649 73513.53125

3.1405 73443.66406

3.13454 73483.40625

3.1286 73644.71875

3.12269 73934.9375

3.1168 74354.00781

3.11093 74892.89844

3.10508 75529.38281

3.09925 76229.77344

3.09345 76945.22656

3.08767 77614.50781

3.08191 78166.33594

3.07618 78523.22656

3.07046 78607.36719

3.06477 78345.64063

3.05909 77675.59375

3.05344 76552.35156

3.04781 74951.97656

3.0422 72870.59375

3.03662 70331.23438

3.03105 67377.13281

3.0255 64069.28516

3.01998 60486.94531

3.01447 56713.59375

3.00899 52842.90625

3.00352 48961.83594

2.99808 45156.04688

2.99265 41497.54688

2.98725 38048.77344

2.98186 34854.71484

2.9765 31949.08398

2.97115 29346.73242

2.96583 27053.16406

2.96052 25059.89063

2.95523 23350.06055

2.94996 21901.12305

2.94472 20684.61914

2.93948 19671.40039

2.93427 18831.74023

2.92908 18137.01563

2.92391 17561.13086

2.91875 17080.80078

2.91362 16676.2168

2.9085 16330.32227

2.9034 16029.45313

2.89831 15762.77734

2.89325 15521.56152

2.88821 15299.38477

2.88318 15091.20898

2.87817 14893.53711

2.87317 14703.64746

2.8682 14519.72363

2.86324 14340.42481

2.8583 14164.88965

2.85338 13992.52148

2.84847 13822.85938

2.84359 13655.73926

2.83871 13490.86426

2.83386 13328.16699

2.82902 13167.64063

2.8242 13009.13769

2.8194 12852.65527

2.81461 12698.16211

2.80984 12545.68066

2.80509 12395.09668

2.80035 12246.43164];

xdata = data(:,1);

ydata = data(:,2);

% Initial guess for parameters:

Eg0 = 3.1;

Eb0 = 0.05;

gama0 = 0.05;

b0 = 0.5;

theta0 = [Eg0; Eb0; gama0; b0];

% lsqcurvefit is in the optimization toolbox.

% fit, from the curve fitting toolbox may be an alternative

theta = lsqcurvefit (@Kdvec, theta0, xdata, ydata);

% check fit;

figure;

plot(xdata, ydata, 'or')

hold on

xplot = linspace(min(xdata),max(xdata));

plot(xplot, Kdvec(theta, xplot))

hold off

end

function Kvec = Kdvec(theta, xdata)

% Vector of Kd for every xdata;

Kvec = zeros(size(xdata));

for i = 1:length(xdata)

Kvec(i) = Kd(theta, xdata(i));

end

end

function K = Kd(theta,rm)

% calculates integral Assumes integrand (1000) is negligible

A=4000;

Eg = theta(1);

Eb = theta(2);

gama = theta(3);

b = theta(4);

integrand = @(r) sech((rm.-r)/gama).* ...

(1/(1-exp(-2*pi*sqrt(Eb/(r-Eg))))).* ...

(1/(1-b*(r-Eg)));

K = A*(1/rm.)*(2*Eb*sech((rm.-Eg+Eb)/gama)+ ...

2*Eb*(1/8)*sech((rm.-Eg+Eb/4)/gama)+ ...

2*Eb*(1/27)*sech((rm.-Eg+Eb/9)/gama)+ ...

integral(integrand,Eg,20);

end

Answer 1

[Torsten]

integrand = @(r) sech((rm-r)/gama).* ...
    (1./(1-exp(-2*pi*sqrt(Eb./(r-Eg))))).* ...
    (1./(1-b*(r-Eg)));
K = A*(1/rm)*(2*Eb*sech((rm-Eg+Eb)/gama)+ ...
    2*Eb*(1/8)*sech((rm-Eg+Eb/4)/gama)+ ...
    2*Eb*(1/27)*sech((rm-Eg+Eb/9)/gama)+ ...
    integral(integrand,Eg,20));

instead of

integrand = @(r) sech((rm.-r)/gama).* ...
    (1/(1-exp(-2*pi*sqrt(Eb/(r-Eg))))).* ...
    (1/(1-b*(r-Eg)));
K = A*(1/rm.)*(2*Eb*sech((rm.-Eg+Eb)/gama)+ ...
    2*Eb*(1/8)*sech((rm.-Eg+Eb/4)/gama)+ ...
    2*Eb*(1/27)*sech((rm.-Eg+Eb/9)/gama)+ ...
    integral(integrand,Eg,20);

[Chaoping]

Thank you, Torsten. It does work, though the fitting was not successful.