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

4 次查看（过去 30 天）
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori 2023-12-27
评论： Alex Sha 2023-12-30
Ran in:
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 个评论
显示 3更早的评论
Sam Chak
Sam Chak 2023-12-29
@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:
https://en.wikipedia.org/wiki/Stopping_power_(particle_radiation)
https://en.wikipedia.org/wiki/Bethe_formula

请先登录，再进行评论。

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

采纳的回答

Alex Sha
Alex Sha 2023-12-28
移动：Sam Chak 2023-12-28
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 个评论
显示 4更早的评论
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori 2023-12-29
Very thanks Dear Alex Sha, if you can provide me code or m.file of the fitting
Alex Sha
Alex Sha 2023-12-30
@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');

请先登录，再进行评论。

更多回答（2 个）

Alex Sha
Alex Sha 2023-12-28
移动：Dyuman Joshi 2023-12-28
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 个评论
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori 2023-12-28
移动：Dyuman Joshi 2023-12-28
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 2023-12-28
移动：Dyuman Joshi 2023-12-28
Ran in:
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.

请先登录，再进行评论。

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023-12-27
Ran in:
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)

类别

在 Help Center 和 File Exchange 中查找有关 Interpolation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by