Compute sensitivity coefficient from numerical gradient

2 vues (au cours des 30 derniers jours)
mz123
mz123 le 21 Nov 2023
Modifié(e) : Matt J le 22 Nov 2023
Hello everybody. I have to rewrite a very basic R script in matlab. The script performs calculation of the sensitivity coefficients, i.e., the partial derivatives, for an example of law of propagation of uncertainty.
Here the R code:
# measurement function
mass2.x <- function(x) {
m.std = x[1]; dm.d = x[2];
diff = x[3]; dm.c = x[4]; dm.B = x[5]
m.std + dm.d + diff + dm.c + dm.B
}
require(numDeriv)
m.std = 10000.005; dm.d = 0.0; diff = mean(c(0.01,0.03,0.02))
dm.c = 0.0; dm.B = 0.0;
sens = grad(func=mass2.x,x=c(m.std,dm.d,diff,dm.c,dm.B))
m.std.u = 0.0225
dm.d.u = 0.015/sqrt(3); dm.c.u = 0.010/sqrt(3)
diff.u = 0.025/sqrt(3); dm.B.u = 0.010/sqrt(3)
m.x = mass2.x(c(m.std,dm.d,diff,dm.c,dm.B))
m.x.unc = sqrt(sum(sens^2*c(m.std.u,dm.d.u,diff.u,dm.c.u,dm.B.u)^2))
basically the measurement model is like
and from the example if I were to calculate the sensitivity coefficient I need to obtain
>> sens = [1 1 1 1 1]
since the model is linear in all the parameters.
Now I tried this approach but I'm not obtaining the same results for sens:
m_std = 10000.005;
dm_d = 0.0;
diff = mean([0.01,0.03,0.02]);
dm_c = 0.0;
dm_B = 0.0;
F = [m_std, dm_d, diff, dm_c, dm_B]
sens = gradient(F)
m_std_u = 0.0225;
dm_d_u = 0.015/sqrt(3);
dm_c_u = 0.010/sqrt(3);
diff_u = 0.025/sqrt(3);
dm_B_u = 0.010/sqrt(3);
m_x = sum(F)
m_x_unc = sqrt(sum(sens.^2 .* [m_std_u,dm_d_u,diff_u,dm_c_u,dm_B_u].^2))
Surely I'm missing something and I do not know if I'm using the function correctly.
Thank you for your help.

Réponses (1)

Matt J
Matt J le 22 Nov 2023
Modifié(e) : Matt J le 22 Nov 2023
I have no background in R. However, if you are trying to take the numerical gradient of a function, there is no native Matlab function that will do that. However this FEX download will,
func=@(x) sum(x);
sens = gradest(func,[0,0,0,0,0])
sens = 1×5
1 1 1 1 1

Catégories

En savoir plus sur Physics dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by