can help me to found empirical equation for data L1 vs T1

6 views (last 30 days)
L1 = [0.784 0.804 0.8845 0.9707 1.062 1.156 1.256 1.36 1.469 1.579 1.696 1.815 1.937 2.063 2.191 2.323 2.457 2.595 2.735 2.876 3.022 3.167 3.316 3.618 3.773 3.927 4.085 4.242 4.405 4.566 4.726 4.891 5.057 5.224 5.236 5.255 5.274 5.242 5.209 5.19 4.634 4.254 3.337 1.257 0.2902 0.11 0.09331 0.02474 0.001703 0.00103 0 0 0];
T1 = [0.9379 1.002 1.25 1.50 1.751 2 2.251 2.502 2.75 3 3.252 3.50 3.752 4.001 4.25 4.501 4.751 5.001 5.253 5.50 5.754 6.002 6.253 6.752 7.003 7.252 7.503 7.75 8.004 8.254 8.5 8.751 9.004 9.256 9.276 9.309 9.375 9.414 9.434 9.44 9.50 9.52 9.546 9.58 9.6 9.609 9.61 9.621 9.643 9.647 10 11 12];
plot(T1, L1, '.-'), grid on
xlabel T_1, ylabel L_1
  5 Comments
Sam Chak
Sam Chak on 29 Dec 2023
@Mushtaq Al-Jubbori, Thanks for the information. Apparently, the Bragg curve describes the stopping power acting on charged particles over the total path length traversed by them, and it is often represented by the Bethe–Bloch formula. It appears that @Alex Sha's formula is more compact. Do you have plans to publish Alex Sha's formula?
References:

Sign in to comment.

Accepted Answer

Alex Sha
Alex Sha on 28 Dec 2023
Moved: Sam Chak on 28 Dec 2023
How about the one below, much simple and works well:
Sum Squared Error (SSE): 0.348116593172766
Root of Mean Square Error (RMSE): 0.0810446642724449
Correlation Coef. (R): 0.999051954884231
Parameter Best Estimate
--------- -------------
p1 1.93075962947385
p2 1.42899205819801
p3 6.09575438791542
p4 -3.221136445554E-15
p5 -4.32095917648694E-12
p6 13.1725280684151
  6 Comments
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori on 29 Dec 2023
Very thanks Dear Alex Sha, if you can provide me code or m.file of the fitting
Alex Sha
Alex Sha on 30 Dec 2023
@Sam Chak: not done in Origin Pro library, but in 1stOpt, there are over ten-thousand functions in its formula library.
@Mushtaq Al-Jubbori: once known the type of fitting function and corresponded data, it is easy to write the fitting code for Matlab like the below, however, the problem now is how to guess the start-values so ensure the convergence, even employ global optimization algorithm like Matlab GA, it is still very difficult to get the right results like ones I provided above. So the most important thing is the tools or algorithms which ensure the computation results are the best solution.
x = [0.9379,1.002,1.25,1.5,1.751,2,2.251,2.502,2.75,3,...
3.252,3.5,3.752,4.001,4.25,4.501,4.751,5.001,5.253,5.5,...
5.754,6.002,6.253,6.752,7.003,7.252,7.503,7.75,8.004,8.254,...
8.5,8.751,9.004,9.256,9.276,9.309,9.375,9.414,9.434,9.44,...
9.5,9.52,9.546,9.58,9.6,9.609,9.61,9.621,9.643,9.647,...
10,11,12];
y = [0.784,0.804,0.8845,0.9707,1.062,1.156,1.256,1.36,1.469,1.579,...
1.696,1.815,1.937,2.063,2.191,2.323,2.457,2.595,2.735,2.876,...
3.022,3.167,3.316,3.618,3.773,3.927,4.085,4.242,4.405,4.566,...
4.726,4.891,5.057,5.224,5.236,5.255,5.274,5.242,5.209,5.19,...
4.634,4.254,3.337,1.257,0.2902,0.11,0.09331,0.02474,0.001703,0.00103,...
0,0,0];
Fit_Fun = @(p,x) (exp(p(1).*x+p(2).*x.^2))./(p(3)-p(4).*exp(0-p(5).*x.^p(6)));
InitialGuess = ones(6,1); %[1,1,1,1,1,1];
p = lsqcurvefit(Fit_Fun, InitialGuess, x, y);
plot(x, y,'pg');
hold on;
plot(x, Fit_Fun(p, x), '-r');
hold off;
grid;
xlabel('x');
ylabel('y');

Sign in to comment.

More Answers (2)

Alex Sha
Alex Sha on 28 Dec 2023
Moved: Dyuman Joshi on 28 Dec 2023
It is not easy to find a function to describe the curve. Refer to the fitting function and result below:
Sum Squared Error (SSE): 0.0168127809048645
Root of Mean Square Error (RMSE): 0.0178107349995405
Correlation Coef. (R): 0.999951618092273
R-Square: 0.999903238525356
Parameter Best Estimate
--------- -------------
p1 -113.045620847505
p2 -1.98381297332703
p3 0.199823219199342
p4 -37.8432173816184
p5 2.28764952811147
p6 1.80019545749072
p7 13.5735683684756
p8 3.19896435789107
p9 1082.4843757759
p10 21.6922662996556
p11 -1.93061546894078
p12 360.858731186603
p13 -8.53668056447937
  2 Comments
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori on 28 Dec 2023
Moved: Dyuman Joshi on 28 Dec 2023
No good, I want one equation for the total curve My be the equation type of tanh, 1/cos or exp
with litle coffetions
John D'Errico
John D'Errico on 28 Dec 2023
Moved: Dyuman Joshi on 28 Dec 2023
Sorry, but just wanting something to be found, or even to exist is not sufficient. (If it was, then I'd have found a nice shiny Lamborghini in my garage wrapped up as a Christmas present.) As you can see, @Alex Sha gave you a shot using a sum of exponentials, but that fails miserably on the rght end, where it starts oscillating. There is no magic way to know what functional form MIGHT fit that set of points. And the very sharp transitions, then changing shape into smooth curves, or well-behaved sections means that any simple functions you choose (certainly NOT polynomials) will fail.
This is why tools like splines (and piecewise polynomials in general) were invented, to allow you to create curves that will represent such things. However, splines are not functions where you can write down the form either. So we can do this:
L1 = [0.784 0.804 0.8845 0.9707 1.062 1.156 1.256 1.36 1.469 1.579 1.696 1.815 1.937 2.063 2.191 2.323 2.457 2.595 2.735 2.876 3.022 3.167 3.316 3.618 3.773 3.927 4.085 4.242 4.405 4.566 4.726 4.891 5.057 5.224 5.236 5.255 5.274 5.242 5.209 5.19 4.634 4.254 3.337 1.257 0.2902 0.11 0.09331 0.02474 0.001703 0.00103 0 0 0];
T1 = [0.9379 1.002 1.25 1.50 1.751 2 2.251 2.502 2.75 3 3.252 3.50 3.752 4.001 4.25 4.501 4.751 5.001 5.253 5.50 5.754 6.002 6.253 6.752 7.003 7.252 7.503 7.75 8.004 8.254 8.5 8.751 9.004 9.256 9.276 9.309 9.375 9.414 9.434 9.44 9.50 9.52 9.546 9.58 9.6 9.609 9.61 9.621 9.643 9.647 10 11 12];
fun = pchip(T1,L1);
fnplt(fun)
hold on
plot(T1,L1,'ro')
grid on
As you can see, it very nicely fits the data. You can evaluate that function at any point, and get a nice smooth looking result. However, there is no simple function you can write down as a result.
fnval(fun,3.75)
ans = 1.9360
The function encapsulated inside fun is actually a set of many tiny segments of cubic polynomials. And that is likely to be the only way to turn that set of points into a function you can then use.

Sign in to comment.


Sulaymon Eshkabilov
Sulaymon Eshkabilov on 27 Dec 2023
Just looking at your data, a first part of it looks like a nice polynomial fit (until it it drops). Here is one easy fit model with poly2, e.g.:
L1 = [0.784 0.804 0.8845 0.9707 1.062 1.156 1.256 1.36 1.469 1.579 1.696 1.815 1.937 2.063 2.191 2.323 2.457 2.595 2.735 2.876 3.022 3.167 3.316 3.618 3.773 3.927 4.085 4.242 4.405 4.566 4.726 4.891 5.057 5.224 5.236 5.255 5.274 5.242 5.209 5.19 4.634 4.254 3.337 1.257 0.2902 0.11 0.09331 0.02474 0.001703 0.00103 0 0 0];
T1 = [0.9379 1.002 1.25 1.50 1.751 2 2.251 2.502 2.75 3 3.252 3.50 3.752 4.001 4.25 4.501 4.751 5.001 5.253 5.50 5.754 6.002 6.253 6.752 7.003 7.252 7.503 7.75 8.004 8.254 8.5 8.751 9.004 9.256 9.276 9.309 9.375 9.414 9.434 9.44 9.50 9.52 9.546 9.58 9.6 9.609 9.61 9.621 9.643 9.647 10 11 12];
plot(T1, L1, '.-'), grid on
xlabel T_1, ylabel L_1
hold on
%% Partial fit because the other part causes the bad fit
IDX = T1<9.5;
T1_S =T1(IDX);
L1_S = L1(IDX);
% Fit Model:
Model = fittype('poly2');
FModel = fit(T1_S', L1_S', Model)
FModel =
Linear model Poly2: FModel(x) = p1*x^2 + p2*x + p3 Coefficients (with 95% confidence bounds): p1 = 0.01782 (0.01579, 0.01986) p2 = 0.3546 (0.3318, 0.3773) p3 = 0.3858 (0.3325, 0.439)
plot(FModel, T1_S, L1_S)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!