Risoluzione di un'equazione differenziale
Metodo di Eulero.

Nell'esempio seguente si considera un'equazione differenziale del primo ordine, nella quale, cioè, compare la funzione e la sua derivata prima.

Un'equazione differenziale si presenta nella forma:

$\tau ·$ \frac{\mathrm{\Delta y(t)}}{\mathrm{\Delta t}}$+\; y(t)\; =$v_i(t)

che, ricondotta nella seguente forma:

$\frac{\mathrm{\Delta y(t)}}{\mathrm{\Delta t}}$ = $\frac{\mathrm{<i>v</i>_i(t)\; -\; y(t)}}{\tau }$

permette di interpretare il problema posto da un'equazione differenziale nei seguenti termini:

determinare il grafico della funzione y(t) conoscendo la tangente alla funzione.
 

La soluzione di un'equazione differenziale è una famiglia di curve parallele. Affinchè la soluzione sia unica bisogna specificare un punto attraverso cui passa la funzione (in tal modo se ne seleziona una). Il punto scelto viene denominato condizione iniziale, e rappresenta il valore della funzione nel punto t=0.

Il metodo di Eulero, quindi, consiste nel calcolare la tangente alla funzione nel punto t=0:

$\frac{\mathrm{\Delta y(0)}}{\mathrm{\Delta t}}$ = $\frac{\mathrm{<i>v</i>_i(0)\; -\; y(0)}}{\tau }$

Moltiplicando entrambi i membri per Δt si ottiene:

Δy(0) = $\frac{\mathrm{<i>v</i>_i(0)\; -\; y(0)}}{\tau }$ ·Δt

che rappresenta, se Δt è molto piccolo rispetto a τ, l'incremento della funzione nel passare dal tempo t=0 al tempo t=0+Δt.

La figura mostra che se si conosce il valore della funzione nel punto iniziale, ad esempio y=0 per t=0, e se si conosce l'espressione della tangente, si può calcolare il coefficiente angolare m della tangente alla funzione in tale punto iniziale, che è dato dal rapporto tra il termine noto meno il valore iniziale della funzione diviso per il coefficiente della derivata: m = $\frac{\mathrm{<i>v</i><sub>i</sub>(0)\; -\; y(0)}}{\tau }$

Quindi, con il valore noto della funzione nel punto t=0, e con l'incremento che la funzione subisce nell'intervallo Δt, si calcola il valore della funzione nel punto t=0+Δt:

y(0+Δt) = y(0) + Δy

Si è così calcolato un altro punto, oltre a quello fissato dalla condizione iniziale, attraverso cui passa la funzione

Dalla figura è anche evidente che il nuovo valore della funzione si viene a trovare sulla retta tangente, mentre il vero valore (P) si trova sulla curva della funzione. Tale errore è tanto più trascurabile quanto più Δt è prossimo a zero.

Programma

Si prepari una funzione che rappresenti il termine noto dell'equazione differenziale. Nella pratica sono di interesse la funzione "Onda quadra", simile ad un segnale digitale, e la funzione "sinusoide", che in base al teorema di Fourier consente di ricostruire una funzione qualsiasi.

Onda Quadra

double vi(double t, double semiperiodo) {
  if (((int)(t/semiperiodo))%2) return 10.0; else return -10.0;
}

Questa funzione riceve, come parametri, il tempo t e la durata di un bit, pari alla metà di un periodo, il quoziente tra questi due numeri, calcolato ad intervalli di tempo minori del semiperiodo, rappresenta la numerazione dei bit. Si decide di assegnare il livello alto ai bit in posizione dispari e il livello basso ai bit in posizione dispari.

Per verificare che questa funzione descrive un'onda quadra si calcoli una tabella di valori in un fissato numero di periodi.

int main() {
  FILE *pf;
  pf = fopen("Eulero.txt", "w");
  
  double Quadra, durataBit, t, dt;
  durataBit = 10.0;
  t = 0.0;
  dt = durataBit/20.0;
  do {
    Quadra = (vi(t, durataBit));
    fprintf(pf, "%f\t%f\n",t, Quadra);
  } while((t += dt) < 10.0*durataBit);
  fclose(pf);
  system("PAUSE");
  return 0;
}

Aprire il file Eulero.txt e osservare che si alternano gruppi di 20 valori 10.0 e gruppi di 20 valori -10.0.

Grafico della funzione

Nel file "Eulero.txt", tramite il menu Modifica, "selezionare tutto". Ancora nel menu "modifica" scegliere il comando "Sostituisci" e sostituire il punto di separazione dei decimali, con la virgola.

Ripetere la selezione di tutto il contenuto del file di testo e copiarlo. Aprire Excel e incollarlo. In Excel, inserire il grafico di questi dati. Si vede che rappresentano un'onda quadra.

Nota: in alcune installazioni l'utente può avere modificato le impostazioni internazionali, usando la virgola come separatore delle migliaia. In questo caso l'operazione di sostituzione del punto con la virgola, nel file di testo, deve essere ripristinata.

la funzione Sinusoide

double vi(double t, double semiperiodo) {
  double f = 1/(2*semiperiodo);
  return 10*sin(6.28*f*t);
}

L'intestazione della funzione corrisponde a quella della funzione onda quadra.

La funzione usa il valore del semiperiodo per calcolare la frequenza della funzione. f = 1/(2*semiperiodo)

Si calcola l'angolo descritto al tempo t, cioè: angolo = 2·π·f·t

Si calcola il seno dell'angolo e lo si moltiplica per l'ampiezza del segnale (10).

Algoritmo di Eulero

Utilizzando, come termine noto, prima l'onda quadra poi la sinusoide, si applica l'algoritmo di Eulero all'equazione differenziale seguente:

$\tau ·$ \frac{\mathrm{\Delta y(t)}}{\mathrm{\Delta t}}$+\; y(t)\; =$v_i(t)

int main(int argc, char *argv[]) {
double t, dt, vc, periodo, tau, deltaVc;
  FILE *pf;
  pf = fopen("Eulero.txt", "w");

Nel main vengono dichiarate le variabili e viene creato il file.

  tau = 10.0;
  vc = 0.0;
  t = dt = tau/20.0;
  periodo=2.0*tau; /* modificare il periodo 5.0*tau, tau/5.0d ecc... */

I valori iniziali assegnati alle variabili hanno il seguente significato:

• il coefficiente della derivata ha il valore 10. La scelta è arbitraria, infatti, si noterà che i risultati ottenuti sono espressi in rapporto a tale valore e quindi il grafico della soluzione sarà una riproduzione in scala della soluzione della stessa equazione differenziale con un diverso valore del coefficiente della derivata.
• il valore della funzione al tempo t=0 vale 0,
• l'intervallo dt si sceglie molto piccolo rispetto al coefficiente τ della derivata prima. (Si osservi che τ ha le dimensioni di un tempo)
• Ci si prepara a calcolare il valore della funzione al tempo t = dt.
• Si sceglie il periodo della funzione di una grandezza paragonabile a τ. Si suggerisce di osservare la soluzione per valori del periodo inferiori a τ (1/2, 1/5, ...) e per valori di τ superiori a τ (2, 5, 10, ...).

  do {
    deltaVc = (vi(t, periodo) - vc)/tau*dt;
    vc += deltaVc;

Inizia l'iterazione per il calcolo dei punti attraverso cui passa la soluzione dell'equazione differenziale.

La prima istruzione di ogni ciclo calcola l'incremento che la funzione subisce nell'intervallo Δt:

ΔVc(t) = $\frac{\mathrm{<i>v</i>_i(t)\; -\; Vc(t)}}{\tau }$ ·Δt

La seconda istruzione calcola il nuovo valore della soluzione, assumendo che se la funzione al tempo t ha il valore vc(t), allora al termine dell'intervallo Δt, in cui è variata della quantià ΔV, ha il valore: Vc(t+Δt) = Vc(t) + ΔV

    fprintf(pf, "%f\t%f\n",t,vc);
  } while((t+=dt) < 5.0*periodo);

I punti calcolati ad ogni ciclo, cioè le coppie (t, vc), vengono registrati nel file di testo, e si passa ad una nuova iterazione incrementando il tempo. Si osserva il comportamento della funzione in 5 periodi della funzione nota.

  fclose(pf);
  system("PAUSE");
  return EXIT_SUCCESS;
}

Problemi

Salvare, insieme al tempo e al valore vc, anche il valore del termine noto vi, separandolo con una tabulazione dagli altri due valori.

Importare i dati in Excel e disegnare il grafico della soluzione. Ripetere questa operazione per vari valori del periodo del termine noto.

Con riferimento alla sinusoide, ipotizzare che 1/τ sia il valore della banda passante e verificare il fenomeno dell'attenuazione, calcolando il rapporto vc(max)/vi(max) per diversi valori della frequenza della funzione nota.