Queste note sono basate sugli appunti fatti con Gianmarco Manzini negli anni 1995-2005

Zeri di funzione (esempi)

Appunti di calcolo numerico

Autore/Autrice
Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Alcuni esempi di problemi che si risolvono trovando le radici di f(x)=0.

Solutori in R

Metodo di bisezione per trovare la radice

bisection_method <- function(func, lower, upper, tol = 1e-6, max_iter = 100) {
  for (i in 1:max_iter) {
    cat(sprintf("Intervallo: lower = %g, upper = %g\n", lower, upper))
    midpoint <- (lower + upper) / 2
    if (func(midpoint) == 0 || (upper - lower) / 2 < tol) {
      return(midpoint)  # Restituisce il punto medio se è la radice
    }
    if (sign(func(midpoint)) == sign(func(lower))) {
      lower <- midpoint  # Aggiorna l'intervallo inferiore
    } else {
      upper <- midpoint  # Aggiorna l'intervallo superiore
    }
  }
  return(NULL)  # Se non si trova la radice
}

Metodo delle secanti per trovare la radice

secant_method <- function(func, x0, x1, tol = 1e-6, max_iter = 100) {
  for (i in 1:max_iter) {
    if (abs(x1 - x0) < tol) {
      return(x1)  # Restituisce la soluzione se la differenza è inferiore alla tolleranza
    }
    x2 <- x1 - func(x1) * (x1 - x0) / (func(x1) - func(x0))  # Calcola il nuovo punto
    cat(sprintf("Iterazione %d: x = %g\n", i, x2))
    x0 <- x1
    x1 <- x2
  }
  warning("La radice non è stata trovata dopo il numero massimo di iterazioni.")
  return(NULL)  # Se non si trova la radice
}

Metodo di Newton-Raphson per trovare la radice

newton_method <- function(func, func_prime, initial_guess, tol = 1e-6, max_iter = 100) {
  x <- initial_guess
  for (i in 1:max_iter) {
    x_new <- x - func(x) / func_prime(x)  # Calcola il nuovo punto
    cat(sprintf("Iterazione %d: x = %g\n", i, x_new))
    if (abs(x_new - x) < tol) {
      return(x_new)  # Restituisce la soluzione se la differenza è inferiore alla tolleranza
    }
    x <- x_new
  }
  warning("La radice non è stata trovata dopo il numero massimo di iterazioni.")
  return(NULL)  # Se non si trova la radice
}

Calcolo della radice cubica

Il calcolo delle radici cubiche è essenziale in ingegneria e in scienze applicate, poiché molte grandezze fisiche e ingegneristiche, come il volume e la densità, richiedono la determinazione di dimensioni lineari a partire da quantità cubiche. In questo documento, analizzeremo il problema di trovare la radice cubica di un numero positivo a risolvendo l’equazione:

f(x) = x^3 - a = 0

Dove a rappresenta il volume di un cubo, e vogliamo determinare la lunghezza del lato x.

Problema

Consideriamo il problema di trovare la radice cubica di un numero a = 27. L’equazione che dobbiamo risolvere è:

f(x) = x^3 - 27 = 0

Stiamo cercando il valore di x tale che x^3 = 27.

Grafico della funzione

Possiamo visualizzare il grafico della funzione per identificare la posizione della radice:

# Grafico della funzione
library(ggplot2)

# Creare un vettore di valori x
x_values <- seq(-5, 5, by = 0.1)
f_values <- x_values^3 - 27

# Creazione del grafico
ggplot(data = data.frame(x = x_values, f = f_values), aes(x = x, y = f)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 3, linetype = "dashed", color = "green") +
  labs(title = "Grafico della funzione $f(x) = x^3 - 27$",
       x = "$x$",
       y = "$f(x)$") +
  theme_minimal()

Soluzione in R

Definizione della funzione per il calcolo della radice cubica

radice_cubica <- function(x) {
  return(x*x*x-27)
}

Metodo di Bisezione

Trova la radice cubica di 27 usando il metodo di bisezione

sol <- bisection_method(radice_cubica, 0, 27)
Intervallo: lower = 0, upper = 27
Intervallo: lower = 0, upper = 13.5
Intervallo: lower = 0, upper = 6.75
Intervallo: lower = 0, upper = 3.375
Intervallo: lower = 1.6875, upper = 3.375
Intervallo: lower = 2.53125, upper = 3.375
Intervallo: lower = 2.95312, upper = 3.375
Intervallo: lower = 2.95312, upper = 3.16406
Intervallo: lower = 2.95312, upper = 3.05859
Intervallo: lower = 2.95312, upper = 3.00586
Intervallo: lower = 2.97949, upper = 3.00586
Intervallo: lower = 2.99268, upper = 3.00586
Intervallo: lower = 2.99927, upper = 3.00586
Intervallo: lower = 2.99927, upper = 3.00256
Intervallo: lower = 2.99927, upper = 3.00092
Intervallo: lower = 2.99927, upper = 3.00009
Intervallo: lower = 2.99968, upper = 3.00009
Intervallo: lower = 2.99989, upper = 3.00009
Intervallo: lower = 2.99999, upper = 3.00009
Intervallo: lower = 2.99999, upper = 3.00004
Intervallo: lower = 2.99999, upper = 3.00001
Intervallo: lower = 2.99999, upper = 3
Intervallo: lower = 2.99999, upper = 3
Intervallo: lower = 3, upper = 3
Intervallo: lower = 3, upper = 3
print(sol)
[1] 3.000001

Metodo delle secanti

Trova la radice cubica di 27 usando il metodo delle secanti

sol <- secant_method(radice_cubica, x0=0, x1=27, max_iter=40)
Iterazione 1: x = 0.037037
Iterazione 2: x = 0.0740232
Iterazione 3: x = 2814.65
Iterazione 4: x = 0.0740266
Iterazione 5: x = 0.07403
Iterazione 6: x = 1642.33
Iterazione 7: x = 0.07404
Iterazione 8: x = 0.07405
Iterazione 9: x = 1641.59
Iterazione 10: x = 0.0740601
Iterazione 11: x = 0.0740701
Iterazione 12: x = 1640.7
Iterazione 13: x = 0.0740801
Iterazione 14: x = 0.0740901
Iterazione 15: x = 1639.81
Iterazione 16: x = 0.0741002
Iterazione 17: x = 0.0741102
Iterazione 18: x = 1638.92
Iterazione 19: x = 0.0741203
Iterazione 20: x = 0.0741303
Iterazione 21: x = 1638.03
Iterazione 22: x = 0.0741404
Iterazione 23: x = 0.0741504
Iterazione 24: x = 1637.14
Iterazione 25: x = 0.0741605
Iterazione 26: x = 0.0741706
Iterazione 27: x = 1636.26
Iterazione 28: x = 0.0741807
Iterazione 29: x = 0.0741908
Iterazione 30: x = 1635.37
Iterazione 31: x = 0.0742008
Iterazione 32: x = 0.0742109
Iterazione 33: x = 1634.48
Iterazione 34: x = 0.074221
Iterazione 35: x = 0.0742312
Iterazione 36: x = 1633.59
Iterazione 37: x = 0.0742413
Iterazione 38: x = 0.0742514
Iterazione 39: x = 1632.7
Iterazione 40: x = 0.0742615
Warning in secant_method(radice_cubica, x0 = 0, x1 = 27, max_iter = 40): La
radice non è stata trovata dopo il numero massimo di iterazioni.
print(sol)
NULL

non funziona proviamo con nuovi punti di partenza

sol <- secant_method(radice_cubica, x0=27/2, x1=27)
Iterazione 1: x = 11.5926
Iterazione 2: x = 10.2912
Iterazione 3: x = 7.33534
Iterazione 4: x = 5.77205
Iterazione 5: x = 4.49521
Iterazione 6: x = 3.69196
Iterazione 7: x = 3.2295
Iterazione 8: x = 3.04379
Iterazione 9: x = 3.00316
Iterazione 10: x = 3.00005
Iterazione 11: x = 3
Iterazione 12: x = 3
print(sol)
[1] 3

Metodo di Newton

Definizione della derivata della funzione

radice_cubica_derivative <- function(x) {
  return(3*x*x)
}

Trova la radice cubica di 27 usando il metodo di Newton

sol <- newton_method(radice_cubica, radice_cubica_derivative, initial_guess = 27)
Iterazione 1: x = 18.0123
Iterazione 2: x = 12.036
Iterazione 3: x = 8.08611
Iterazione 4: x = 5.52838
Iterazione 5: x = 3.98006
Iterazione 6: x = 3.22152
Iterazione 7: x = 3.01488
Iterazione 8: x = 3.00007
Iterazione 9: x = 3
Iterazione 10: x = 3
print(sol)
[1] 3

Equazione di radiazione di Stefan-Boltzmann (Trasferimento di calore)

Problema

La potenza irradiata da un corpo nero è data dalla legge di Stefan-Boltzmann:

P = \sigma A (T^4 - T_{\infty}^4)

Dove: - P = 500 \, \text{W} - A = 2 \, \text{m}^2 - T_{\infty} = 300 \, \text{K} - \sigma = 5.67 \times 10^{-8} \, \text{W/m}^2\text{K}^4

Trova la temperatura del corpo T risolvendo l’equazione:

f(T) = 5.67 \times 10^{-8} \cdot 2 \cdot (T^4 - 300^4) - 500

Soluzione in R

Definizione della funzione f(T):

f <- function(T) {
  P <- 500
  A <- 2
  sigma <- 5.67e-8
  return(sigma * A * (T^4 - 300^4) - P)
}

Metodo di Bisezione

Trova la temperatura usando il metodo di bisezione:

sol <- bisection_method(f, 300, 600)
Intervallo: lower = 300, upper = 600
Intervallo: lower = 300, upper = 450
Intervallo: lower = 300, upper = 375
Intervallo: lower = 300, upper = 337.5
Intervallo: lower = 318.75, upper = 337.5
Intervallo: lower = 328.125, upper = 337.5
Intervallo: lower = 332.812, upper = 337.5
Intervallo: lower = 332.812, upper = 335.156
Intervallo: lower = 333.984, upper = 335.156
Intervallo: lower = 333.984, upper = 334.57
Intervallo: lower = 334.277, upper = 334.57
Intervallo: lower = 334.424, upper = 334.57
Intervallo: lower = 334.424, upper = 334.497
Intervallo: lower = 334.424, upper = 334.46
Intervallo: lower = 334.424, upper = 334.442
Intervallo: lower = 334.424, upper = 334.433
Intervallo: lower = 334.428, upper = 334.433
Intervallo: lower = 334.431, upper = 334.433
Intervallo: lower = 334.431, upper = 334.432
Intervallo: lower = 334.431, upper = 334.432
Intervallo: lower = 334.431, upper = 334.432
Intervallo: lower = 334.431, upper = 334.432
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
Intervallo: lower = 334.431, upper = 334.431
print(sol)
[1] 334.4315

Metodo delle secanti

Trova la soluzione per la temperatura usando il metodo delle secanti:

sol <- secant_method(f, x0 = 500, x1 = 400 )
Iterazione 1: x = 364.523
Iterazione 2: x = 341.536
Iterazione 3: x = 335.308
Iterazione 4: x = 334.459
Iterazione 5: x = 334.432
Iterazione 6: x = 334.431
Iterazione 7: x = 334.431
print(sol)
[1] 334.4315

Metodo di Newton

Derivata della funzione f(T) rispetto a T

f_prime <- function(T) {
  A <- 2
  sigma <- 5.67e-8
  return(4 * sigma * A * T^3)
}

Trova la soluzione per la temperatura usando il metodo di Newton-Raphson

sol <- newton_method(f, f_prime, initial_guess = 500)
Iterazione 1: x = 400.018
Iterazione 2: x = 348.871
Iterazione 3: x = 335.303
Iterazione 4: x = 334.435
Iterazione 5: x = 334.431
Iterazione 6: x = 334.431
print(sol)
[1] 334.4315

Equazione di Colebrook (Flusso di fluidi)

Problema

L’equazione di Colebrook per il coefficiente di attrito f è definita come:

\frac{1}{\sqrt{f}} = -2 \log \left( \frac{\epsilon / D}{3.7} + \frac{2.51}{Re \cdot \sqrt{f}} \right)

Dove:

  • \epsilon / D = 0.005
  • Re = 1.5 \times 10^5

Trova il coefficiente di attrito f.

Soluzione in R

Definizione della funzione di Colebrook

colebrook <- function(f) {
  epsilon_D <- 0.005
  Re <- 1.5e5
  return(1/sqrt(f) + 2 * log10(epsilon_D/3.7 + 2.51/(Re * sqrt(f))))
}

Metodo di Bisezione

Trova il coefficiente di attrito usando il metodo di bisezione

sol <- bisection_method(colebrook, 0.01, 0.1)
Intervallo: lower = 0.01, upper = 0.1
Intervallo: lower = 0.01, upper = 0.055
Intervallo: lower = 0.01, upper = 0.0325
Intervallo: lower = 0.02125, upper = 0.0325
Intervallo: lower = 0.026875, upper = 0.0325
Intervallo: lower = 0.0296875, upper = 0.0325
Intervallo: lower = 0.0296875, upper = 0.0310938
Intervallo: lower = 0.0303906, upper = 0.0310938
Intervallo: lower = 0.0307422, upper = 0.0310938
Intervallo: lower = 0.030918, upper = 0.0310938
Intervallo: lower = 0.030918, upper = 0.0310059
Intervallo: lower = 0.0309619, upper = 0.0310059
Intervallo: lower = 0.0309839, upper = 0.0310059
Intervallo: lower = 0.0309949, upper = 0.0310059
Intervallo: lower = 0.0310004, upper = 0.0310059
Intervallo: lower = 0.0310004, upper = 0.0310031
Intervallo: lower = 0.0310017, upper = 0.0310031
print(sol)
[1] 0.03100243

Metodo delle secanti

Trova il coefficiente di attrito usando il metodo delle secanti

sol <- secant_method(colebrook, x0 = 0.02, x1 = 0.03)
Iterazione 1: x = 0.030725
Iterazione 2: x = 0.0309953
Iterazione 3: x = 0.031002
Iterazione 4: x = 0.031002
print(sol)
[1] 0.03100205

Metodo di Newton

Definizione della derivata della funzione di Colebrook rispetto a f

colebrook_prime <- function(f) {
  epsilon_D <- 0.005
  Re <- 1.5e5
  return(-0.5 * f^(-1.5) - (2 * 2.51 / (Re * (f^(1.5)) * log(10))) / (epsilon_D / 3.7 + 2.51 / (Re * sqrt(f))))
}

Trova il coefficiente di attrito usando il metodo di Newton-Raphson

sol <- newton_method(colebrook, colebrook_prime, initial_guess = 0.02)
Iterazione 1: x = 0.0277966
Iterazione 2: x = 0.0307201
Iterazione 3: x = 0.0309974
Iterazione 4: x = 0.031002
Iterazione 5: x = 0.031002
print(sol)
[1] 0.03100205

Equazione di Van der Waals (Equilibrio chimico)

Problema

L’equazione di Van der Waals è data da:

\left( P + \frac{a}{V^2} \right) (V - b) = RT

Dove:

  • P = 100\, \text{kPa}
  • T = 300\, \text{K}
  • a = 0.9\, \text{Pa} \cdot \text{m}^6 / \text{mol}^2
  • b = 0.1\, \text{m}^3 / \text{mol}

Trova il volume molare V.

Soluzione in R

Definizione della funzione di Van der Waals

vanderwaals <- function(V) {
  P <- 100e3  # kPa to Pa
  R <- 8.314  # J/(mol K)
  T <- 300    # K
  a <- 0.9    # Pa m^6/mol^2
  b <- 0.1    # m^3/mol
  return((P + a/(V^2)) * (V - b) - R * T)
}

Metodo di Bisezione

Trova il volume usando il metodo di bisezione

sol <- bisection_method(vanderwaals, 0.1, 10)
Intervallo: lower = 0.1, upper = 10
Intervallo: lower = 0.1, upper = 5.05
Intervallo: lower = 0.1, upper = 2.575
Intervallo: lower = 0.1, upper = 1.3375
Intervallo: lower = 0.1, upper = 0.71875
Intervallo: lower = 0.1, upper = 0.409375
Intervallo: lower = 0.1, upper = 0.254688
Intervallo: lower = 0.1, upper = 0.177344
Intervallo: lower = 0.1, upper = 0.138672
Intervallo: lower = 0.119336, upper = 0.138672
Intervallo: lower = 0.119336, upper = 0.129004
Intervallo: lower = 0.12417, upper = 0.129004
Intervallo: lower = 0.12417, upper = 0.126587
Intervallo: lower = 0.12417, upper = 0.125378
Intervallo: lower = 0.124774, upper = 0.125378
Intervallo: lower = 0.124774, upper = 0.125076
Intervallo: lower = 0.124925, upper = 0.125076
Intervallo: lower = 0.124925, upper = 0.125001
Intervallo: lower = 0.124925, upper = 0.124963
Intervallo: lower = 0.124925, upper = 0.124944
Intervallo: lower = 0.124925, upper = 0.124935
Intervallo: lower = 0.124925, upper = 0.12493
Intervallo: lower = 0.124928, upper = 0.12493
Intervallo: lower = 0.124928, upper = 0.124929
print(sol)
[1] 0.1249282

Metodo delle secanti

Trova il volume usando il metodo delle secanti

sol <- secant_method(vanderwaals, x0 = 0.5, x1 = 1)
Iterazione 1: x = 0.124923
Iterazione 2: x = 0.124928
Iterazione 3: x = 0.124928
print(sol)
[1] 0.1249276

Metodo di Newton

Derivata della funzione di Van Der Waals

vanderwaals_derivative <- function(V) {
  P <- 100e3  # Pa
  R <- 8.314  # J/(mol K)
  T <- 300    # K
  a <- 0.9    # Pa m^6/mol^2
  b <- 0.1    # m^3/mol

  # Calcolo della derivata
  term1 <- P + a / (V^2)
  term2 <- (2 * a * (V - b)) / (V^3)

  return(term1 - term2)
}

Trova il volume usando il metodo di Newton

sol <- newton_method(vanderwaals, vanderwaals_derivative, initial_guess = 1)
Iterazione 1: x = 0.124928
Iterazione 2: x = 0.124928
print(sol)
[1] 0.1249276

Equazione di Shockley per i diodi

Problema

La corrente che attraversa un diodo è descritta dall’equazione di Shockley: I = I_0 \left( e^{\frac{V}{nV_T}} - 1 \right) Dove: - I = 2 \, \text{mA} (corrente attraverso il diodo) - I_0 = 10^{-12} \, \text{A} (corrente di saturazione inversa) - n = 1 (fattore di idealità del diodo) - V_T = 26 \, \text{mV} (tensione termica)

Trova la tensione V del diodo.

Soluzione in R

# Definizione della funzione di Shockley
shockley <- function(V) {
  I <- 2e-3  # Corrente attraverso il diodo in Ampere
  I0 <- 1e-12  # Corrente di saturazione inversa in Ampere
  n <- 1  # Fattore di idealità del diodo
  VT <- 26e-3  # Tensione termica in Volt
  return(I0 * (exp(V / (n * VT)) - 1) - I)
}

Metodo di Bisezione

Trova la tensione del diodo usando il metodo di bisezione

sol <- bisection_method(shockley, 0, 1)
Intervallo: lower = 0, upper = 1
Intervallo: lower = 0.5, upper = 1
Intervallo: lower = 0.5, upper = 0.75
Intervallo: lower = 0.5, upper = 0.625
Intervallo: lower = 0.5, upper = 0.5625
Intervallo: lower = 0.53125, upper = 0.5625
Intervallo: lower = 0.546875, upper = 0.5625
Intervallo: lower = 0.554688, upper = 0.5625
Intervallo: lower = 0.554688, upper = 0.558594
Intervallo: lower = 0.556641, upper = 0.558594
Intervallo: lower = 0.556641, upper = 0.557617
Intervallo: lower = 0.556641, upper = 0.557129
Intervallo: lower = 0.556641, upper = 0.556885
Intervallo: lower = 0.556763, upper = 0.556885
Intervallo: lower = 0.556824, upper = 0.556885
Intervallo: lower = 0.556824, upper = 0.556854
Intervallo: lower = 0.556824, upper = 0.556839
Intervallo: lower = 0.556824, upper = 0.556831
Intervallo: lower = 0.556824, upper = 0.556828
Intervallo: lower = 0.556826, upper = 0.556828
print(sol)
[1] 0.5568266

Metodo delle secanti

Trova la tensione del diodo usando il metodo delle secanti

sol <- secant_method(shockley, x0 = 0.6, x1 = 0.7)
Iterazione 1: x = 0.598232
Iterazione 2: x = 0.596581
Iterazione 3: x = 0.576856
Iterazione 4: x = 0.567524
Iterazione 5: x = 0.560234
Iterazione 6: x = 0.557467
Iterazione 7: x = 0.556868
Iterazione 8: x = 0.556827
Iterazione 9: x = 0.556827
print(sol)
[1] 0.5568267

Metodo di Newton

Derivata della funzione di Shockley

shockley_derivative <- function(V) {
  I0 <- 1e-12  # Corrente di saturazione inversa in Ampere
  n <- 1       # Fattore di idealità del diodo
  VT <- 26e-3  # Tensione termica in Volt
  exp_term <- exp(V / (n * VT))  # Calcola l'esponenziale
  derivative <- (I0 / (n * VT)) * exp_term  # Calcola la derivata
  return(derivative)
}

Trova la tensione del diodo usando il metodo di Newton

sol <- newton_method(shockley, shockley_derivative, initial_guess = 0.6)
Iterazione 1: x = 0.578941
Iterazione 2: x = 0.564048
Iterazione 3: x = 0.557743
Iterazione 4: x = 0.556843
Iterazione 5: x = 0.556827
Iterazione 6: x = 0.556827
print(sol)
[1] 0.5568267

Bibliografia

  1. F.P. Incropera, D.P. DeWitt, T.L. Bergman, A.S. Lavine.
    Fundamentals of Heat and Mass Transfer.
    John Wiley & Sons. (2006)

  2. F.M. White.
    Fluid Mechanics.
    McGraw-Hill Education. (2011)

  3. J.M. Smith, H.C.Van Ness, M.M. Abbott, M.T. Swihart.
    Introduction to Chemical Engineering Thermodynamics.
    McGraw-Hill Education. (2005)

  4. A.S. Sedra, K.C. Smith, T.C. Carusone, V.Gaudet.
    Microelectronic Circuits.
    Oxford University Press. (2020)