Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- nu = 51
- sig = 0.0011
- K=23490
- mu=3.9
- ### calcul de Emid
- #fonction pour l'intégrale double
- f <- function(t,s){
- a <- (exp(s)-1)/s
- b <- (sig**2*t/6 + sig**2)*exp(-mu*t)*a
- c <- nu*exp(-nu*t)*b
- return(c)
- }
- #fonction pour l'intégrale simple
- g <- function(t){
- a <- nu*exp(-nu*t)
- b <- sig**2*t/3*(1-exp(-mu*t))
- return(a*b)
- }
- emid <- function(h){
- tmin <- h
- tmax <- 100
- smax <- function(x) mu*x
- smin <- 0
- I <- integral2(f,tmin,tmax,smin,smax)$Q
- J <- integrate(g,h,100)$value
- a <- (1 - exp(-nu*h)*(nu/(nu+h)))*exp(nu*h)
- return(a*(I+J))
- }
- #### calcul de E final
- # fonction pour l'intégrale triple
- k <- function(t,s,u){
- a <- (exp(u)-1)/u
- b <- 1/s
- c <- nu*exp(-nu*t)*exp(-mu*t)
- return(a*b*c)
- }
- efin <- function(h){
- tmin <- h
- tmax <- 100
- smax <- function(x) mu*x
- smin <- 0
- umax <- function(x,y) y
- umin <- 0
- I <- integral3(k,tmin,tmax,smin,smax,umin,umax)
- a <- (1-exp(-mu*h)*mu/(mu+h))
- b <- exp(mu*h)/a
- J <- nu**2*(nu*h+0.5)
- e <- emid(h) + (J/K**2)*b*I
- return(e)
- }
- #Graphique
- c <- 0.1 * 1:140
- d <- lapply(c,efin)
- plot(c,d,xlab = 'h en secondes',ylab = 'E(h)',main = 'Courbe de E(h)')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement