Je m'intéresse à la dissolution d'un film d'épaisseur L dans une solution aqueuse. Je recherche le coefficient de diffusion D
J'ai réalisé ce code naïf pour ajuster le paramètre D par rapport à mes données sans succès.
J'ai voulu minimiser l'erreur quadratique.
Code : Tout sélectionner
# Mes données
t <- c(0.0000, 0.0833, 0.1660, 0.2500,
0.3330, 0.4166, 0.5000, 0.6660,
0.8333, 1.0000, 1.5000, 2.0000,
2.5000, 3.0000, 3.5000, 4.0000,
4.5000, 5.0000, 5.5000, 6.0000,
6.5000, 7.0000, 7.5000, 8.0000,
24.0000)
y <- c(0.0000000, 0.1891365, 0.2807108,
0.3522929, 0.3810654, 0.4089780,
0.4273398, 0.4624385, 0.5069652,
0.5271696, 0.5899072, 0.6566372,
0.6942848, 0.7397635, 0.7518605,
0.8004716, 0.8110638, 0.8310532,
0.8723246, 0.8702028, 0.8864763,
0.9042854, 0.9139869, 0.9285407,
1.0000000)
data <- data.frame(t = t, y = y)
# Ma fonction d'ajustement
f <- function(t, D) {
nbtermessomme <- 1000
n <- 0:nbtermessomme
L <- 0.0126
terms <- exp(-D * (2 * n + 1)^2 * pi^2 * t / L^2) / (2 * n + 1)^2
y <- 1 - 8 / pi^2 * sum(terms)
return(y)
}
Dseq = seq(0,1e-4,1e-10)
res=0*Dseq
for (index in 1:length(Dseq))
D_value <- Dseq[index]
result <- sapply(t, function(t) f(t, D_value))
datatemp=cbind(data,result)
res[index]=mean((datatemp$y-datatemp$result)^2)
Dseq[which.min(res)]
Je trouve D=0 or D désigne un paramètre de diffusion et je suis surpris de trouver cette valeur. De plus, y a-t-il une méthode moins naïve pour estimer D ?
Quelqu'un pourrait-il m'aider, s'il vous plaît ?