Teoria 3

Da CdM_unimore.

FORTRAN: Soluzioni di un'equazione di 2°grado (pag.19 P.A.)

Si desidera scrivere un programma in Fortran per la risoluzione di un'equazione di secondo grado, con discussione del discriminante. Data dunque un'equazione di secondo grado nell'incognita x

definito il discriminante di tale funzione come:

E' possibile dunque scrivere il programma per la sua risuoluzione, che nomineremo EQ2

 C    PROGRAM EQ2
      EQUAZIONE DI 2° GRADO a*x**2+b*x+c=0
      EPSI=1.E-5
      print*,'scrivi a,b,c ',A,B,C
      read*, A,B,C

Dove si può osservare come al posto del comando print*,, si potrebbe utilizzare il comando write (*,*), con la sola differenza che mentre il comando print*,, non ha controllo della posizione dell'output sullo schermo, il comando write (*,*) permette di incolonnare. Le lettere A,B,C in dichiarazione implicita indicano l'utilizzo di numeri reali. Il proseguo del codice risulta essere:

      DISCR=B**2-4*A*C
      IF(DISCR.GT.EPSI)GO TO 100
      IF(DISCR.LT.-EPSI) GO TO 200
      IF(ABS(DISCR).LE.EPSI)GO TO 300

Nel maneggiare il discriminante è opportuno non introdurre subito il comando SQRT (radice quadrata),in preludio alla risoluzione dell'equazione, in quanto se il radicando risultasse essere negativo si determinerebbe un segnale di errore. La discussione del segno del determinante e dei casi risolutivi è operata per mezzo della funzione IF. Nella funzione IF le espressioni logiche elementari utilizzabili sono, date due variabili A e B

A.EQ.B vero se A=B

A.NE.B vero se A≠ B

A.GT.B vero se A>B

A.GE.B vero se A≥B

A.LT.B vero se A<B

A.LE.B vero se A≤B

Con questa scrittura, se la funzione discriminante è maggiore di epsi (un numero molto piccolo positivo) il codice successivamente letto è quello presente nella riga etichettata come 100, se invece è minore di - epsi (un valore molto piccolo negativo)si rimanda alla riga di codice 200, mentre se il valore del determinante risulta essere in termini assoluti inferiore a quello di epsi il comando GOTO rinvia alla riga di codice 300. 100, 200, 300 sono quindi delle etichette (o label) che rappresentano determinate zone di un programma. In generale la funzione GOTO ha come conseguenza un "salto", vale a dire, l'esecuzione viene continuata dalla riga di etichetta e non con l'istruzione della riga successiva al GOTO stesso. Se il GOTO non fosse preceduto dall'ipotesi contenuta nell'IF si parlerebbe di un "GOTO incondizionato", poiché il salto viene eseguito sempre. Quando il GOTO è unito con un IF si parla di "GOTO condizionato"; questa istruzione è usata frequentemente per formare loops senza ricorrere all'istruzione DO. In questo caso, se il salto viene eseguito ad una istruzione di una riga precedente, si corre il rischio di creare un loop infinito, cioè un loop senza vie di uscita. Un altra istruzione che può essere utile è il cosiddetto "GOTO calcolato", che ha la seguente forma

GOTO(l1,l2,...,ln),espressione numerica intera

e causa un trasferimento di controllo alla riga etichettata l1 se l'espressione numerica intera valesse 1, alla riga etichettata l2 se l'espressione valesse 2, e così via. Se l'espressione numerica è negativa, zero o maggiore di n, il GOTO calcolato viene semplicemente saltato.

Tornando a considerare il nostro codice per la risoluzione dell'equazione di secondo grado, si può notare inoltre come l'utilizzo di EPSI, valore molto piccolo, è funzionale al fatto che discutere del segno del discriminante, in particolare nel caso di discriminante nullo, rispetto al valore reale 0., potrebbe essere fuorviante. La presenza inevitabile dello "sporco di macchina" potrebbe infatti far corrispondere anche a casi di determinanti nulli i casi di determinante di diverso da zero. L'utilizzo di un epsilon EPSI molto piccolo nella discussione del segno del determinante permette dunque di scongiurare questa eventualità.

      DISCR=B**2-4*A*C
      IF(DISCR.GT.EPSI)GO TO 100
      IF(DISCR.LT.-EPSI) GO TO 200
      IF(ABS(DISCR).LE.EPSI)GO TO 300
      PRINT*,'PROBABILE ERRORE'

Si può inoltre fare in modo, attraverso l'ultima riga inserita, che nel caso le tre condizioni non siano soddisfatte venga visualizzato un messaggio a video di errore, determinando una condizione di "diagnostica interna" al nostro programma. Da notare infine come il terzo comando IF, nella serie di comandi logici, potrebbe essere sostituito da un comando di tipo ELSE, che senza discutere del valore assoluto del discriminante ed escludendo i due casi coperti dai primi due IF, invii direttamente alla riga 300 del codice (quella relativa appunto al caso di determinante nullo). Questa scelta non permetterebbe tuttavia di determinare la "diagnostica interna" così come realizzata di sopra.

Dopo di che è necessario impostare il contenuto delle righe 100, 200, 300. A partire dalla riga 100, corrispondente alla condizione di determinante maggiore di zero, vengono determinate le due soluzioni reali distinte.

  100 continue
      SOL1=(-B+SQRT(DISCR))/(2*A)
      SOL2=(-B-SQRT(DISCR))/(2*A)
      print*, 'Due soluzioni reali'
      print*, 'Sol_1= ',SOL1,' Sol_2= 'SOL2
      stop

La riga 200, corrispondente alla condizione di determinante negativo, informa invece che in questo caso non esistono soluzioni reali.

  200 continue
      print*, 'Nessuna Sol reale'
      stop

La riga 300, relativa alla condizione di determinante nullo, determina infine le due soluzioni reali coincidenti.

  300 continue
      print*, 'La Sol reale è doppia'
      SOL=-B/2/A
      print*, 'Sol= ',SOL
      stop
      end

Codice alternativo

c234567---------------------------------------------------------------!
      program eq2grado2
   	  real epsi, discr, sol1, sol2, a, b, c
   	  
   	  epsi=1.e-6
   	  write(*,*) 'inserire i coefficienti ax^2+bx+c'
   	  read(*,*) a,b,c
   	  discr=b**2-4*a*c
   	  if (discr .gt. epsi) then
	     write(*,*) 'le soluzioni sono reali e distinte'
   		 sol1=(-b+sqrt(discr))/(2*a)
   		 sol2=(-b-sqrt(discr))/(2*a)
   		 write(*,*) 'sol1= ', sol1
   		 write(*,*) 'sol2= ', sol2
   	  elseif (discr .lt. epsi) then
	     write(*,*) 'le soluzioni sono inesistenti'
   	  else
	     write(*,*) 'le soluzioni sono reali e coincidenti'
   		 sol1=(-b+sqrt(discr))/(2*a)
   		 write(*,*) 'sol1= ', sol1
   	  endif
   	  stop 
   	  end

Calcolo sigma-ideale secondo Mohr (es. 4.1.3 pag.19 P.A.)

cerchi di Mohr

Si desidera realizzare un programma in codice Fortran che permetta di calcolare la tensione ideale SIGMAID secondo Mohr (Guest), assegnati i valori delle tensioni principali (SIGMA1, SIGMA2, SIGMA3), il cui segno sia indifferentemente concorde o discorde.

Si chiama il programma in questione PROGRAM MOHR

      print*, 'Batti le tre tensioni reali principali'
      read(*,*) SIGMA1, SIGMA2, SIGMA3
      SIGMAID = AMAX1( ABS( SIGMA1-SIGMA2 ), ABS( SIGMA2-SIGMA3 ), ABS( SIGMA3-SIGMA1 ) )
      print*, 'Sigma_id: ',SIGMAID
      stop
      end

Il ruolo principale all'interno del codice è svolto dalla funzione intrinseca

      AMAX1

la quale permette di determinare il massimo tra i componenti di un vettore.All'interno di Fortran esiste infatti una serie di functions precostituite e che servono per il calcolo delle funzioni matematiche più comuni, come le funzioni trigonometriche, ecc. Per ogni caso, ne esistono diverse versioni a seconda del tipo dell' argomento e della precisione con cui vogliamo conoscere il risultato.Come regola quasi generale si ha che la function che esegue il calcolo in doppia precisione si ottiene anteponendo al nome "normale" la lettera D, per la precisione quadrupla la lettera Q, per avere il risultato con i complessi la lettera C, ecc. In particolare per la nostra funzione possiamo dunque osservare a cosa corrisponderebbero le diverse scritture:

MAX0(X1,.,XN) interi I4 intero I4 massimo tra X1,..,XN;

AMAX1(X1,.,XN) reali R4 reale R4 massimo tra X1,..,XN;

DMAX1(X1,.,XN) reali R8 reali R8 massimo tra X1,..,XN;

Subroutine CLEAR per l'inizializzazione di una matrice (pag.108 P.A.)

Una subroutine è la forma migliore per spezzare il programma in parti più piccole, ognuna dedicata ad un compito specifico. Il programma chiamante accede alla subroutine tramite l' istruzione CALL:

CALL SUNAME(ARG1,...,ARGN)

che causa l' esecuzione del gruppo di istruzioni che iniziano con

SUBROUTINE SUNAME(VAR1,...,VARN)

e finiscono con l' istruzione END. Il controllo ritorna al programma chiamante quando la subroutine termina l' esecuzione con l' istruzione RETURN. Da notare che "SUNAME" è semplicemente un nome usato per chiamare la subroutine e non una variabile. I risultati del calcolo sono passati al programma chiamante attraverso gli argomenti ARG1,...,ARGN della subroutine.

In particolare è nostra intenzione ora realizzare la stesura di una subroutine che permetta l'inizializzazione a zero degli elementi di una matrice. L'inzializzazione degli array o matrici, sebbene sia svolta in automatico in alcuni compilatori, è un'operazione importante in quanto permette di essere sicuri che il valore delle singole posizioni dell'array sia inizialmente pari allo zero macchina evitando quindi errori. Considerare le posizioni di un'array pari allo zero macchina, nonostante esso sia stato dichiarato come 'nuovo' all'inizio del programma, sarebbe infatti sbagliato e alcune posizioni potrebbero conservare i valori assegnate per le stesse celle di memoria all'interno di programmi precedenti. Inizializzare le matrici risulta quindi essere un'operazione fondamentale per evitare errori nei risultati altrimenti difficili da individuare senza questa precauzione.

La matrice da inizializzare è la A, ed i suoi indici sono K1 e K2. Se uno dei due indici vale 1, la matrice diventa un vettore. I valori degli indici dovranno quindi essere necessariamente specificati nel main a cui a è richiamata la subroutine.

SOUBROUTINE CLEAR(K1,K2,A)

K1,K2 sono numeri interi mentre A è una matrice di numeri reali. E' prassi che prima si inseriscano i numeri e poi le matrici. Questi numeri servono e il rispetto di tale ordine nella chiamata del main è necessario.

      SOUBROUTINE CLEAR(K1,K2,A)
 C    pulizia matrice
      DIMENSION A(K1,K2)

Quest'ultimo comando permette di dare una dimensione alla matrice A. Poichè la matrice A è caratterizzata da due indici, per inizializzarla sarà necessario utilizzare due cicli DO:

      SOUBROUTINE CLEAR(K1,K2,A)
 C    pulizia matrice
      DIMENSION A(K1,K2)
      DO10, I10=1,K1
      DO20, I20=1,K2

dove se scrivessi a questo punto

      A(K1,K2)=0

è un errore perchè non pemetterei il corretto scorrimento degli indici. Il modo corretto di scrivere è il seguente:

 C    pulizia matrice
      DIMENSION A(K1,K2)
      DO10, I10=1,K1
         DO20, I20=1,K2
            A(I10,I20)=0
   20    CONTINUE
   10 CONTINUE
      RETURN
      END

Con il return, una volta finita questa subroutine, il programma si riporta al main.

Subroutine di trasposizione di una matrice (pag.109 P.A.)

Così come effettuato nel paragrafo precedente, si intende ora riportare una routine che effettui calcoli di utilità generale, e on calcoli specifici di un programma agli elementi finiti. In particolare si vuole realizzare ora una routine che permetta di effettuare la trasposizione di una matrice. La matrice trasposta di una matrice si ottiene scambiandone le righe con le colonne. La trasposta di una matrice dunque è la matrice il cui generico elemento con indici è l'elemento con indici della matrice originaria. In simboli:

con lo spazio vettoriale delle matrici di dimensione n. In pratica, la matrice trasposta si deve intendere come una matrice in cui le colonne diventano righe e le righe diventano colonne. L'operazione di trasposizione è definita sia su matrici quadrate che rettangolari, e quindi anche su vettori. In particolare, un vettore colonna trasposto è un vettore riga, e viceversa. Una matrice che coincide con la propria trasposta è detta matrice simmetrica, e deve essere una matrice quadrata. Quindi, sebbene in generale date due matrici e di dimensioni opportune si ha che:

l'operatore di trasposizione è una trasformazione lineare, ovvero, dati due scalari ed , vale:

Più in generale, dati N scalari ed N matrici di pari dimensioni, vale:

dove indica una sommatoria.

Indicando nel codice frtran la matrice iniziale con A e la sua trasposta con AT, la subroutine necessaria per ottenere dunque la trasposta di una matrice risulta essere la seguente.

      SUBROUTINE TRSP(K1,K2,A,AT)
 C    trasposizione della matrice A
      DIMENSION A(K1,K2),AT(K2,K1)
      DO 10, I10=1,K1
         DO 20, I20=1,K2
            AT(I20,I10)=A(I10,I20)
   20    CONTINUE
   10 CONTINUE
      RETURN
      END

Dove gli indici della matrice A sono K1 e K2, mentre gli indici della matrice AT sono rispettivamente quindi K2 e K1.

Subroutine matrice D, caso tensione piana (pag.116 P.A.)

Assumendo che un componente meccanico lavori in elasticità lineare, le tensioni all'interno dell'elemento finito triangolare si ottengono dalle deformazioni attraverso la legge di Hooke. Assumendo uno stato di tensione piana, le tensioni nel piano valgono:

dove G è definita come:

In caso di tensione piana la matrice D viene definita quindi nel seguente modo ed esprime in sostanza la legge di Hooke

Anche le tensioni, come le deformazioni, rimangono costanti all'interno dell'elemento finito triangolare, dato che le tensioni si ricavano dalle deformazioni, moltiplicandole per le costanti elastche. La subroutine in Fortran necessaria per compilare la matrice D sarà dunque:

 SUBROUTINE DMAT(YM,PR,D)
 C     LEGGE DI HOOKE
       DIMENSION D(3,3)
       COST1=YM/(1-PR**2)
       D(1,1)=COST1
       D(1,2)=COST1*PR
       D(2,1)=COST1*PR
       D(2,2)=COST1
       D(1,3)=0.
       D(2,3)=0.
       D(3,1)=0.
       D(3,2)=0.
       D(3,3)=YM/(2*(1-PR**2))
       RETURN
       END

Subroutine matrice B (pag.118 P.A.)

La matrice B permette di realizzare la seguente uguaglianza

o in maniera espansa:

Fig4.jpg








Utilizziamo una subroutine di Fortran per scrivere la matrice B

SOUBROUTINE BMAT(XI,YI,XJ,YJ,XK,YK,B)
      DIMENSION B(3,6)
      B(1,1)=YJ-YK
      B(1,2)=0.
      B(1,3)=YK-YI
      B(1,4)=0.
      B(1,5)=YI-YJ
      B(1,6)=0.
      B(2,1)=0.
      B(2,2)=XK-XJ
      B(2,3)=0.
      B(2,4)=XI-XK
      B(2,5)=0.
      B(2,6)=XJ-XI
      B(3,1)=XK-XJ
      B(3,2)=YJ-YK
      B(3,3)=XI-XK
      B(3,4)=YK-YI
      B(3,5)=XJ-XI
      B(3,6)=YI-YJ

A questo punto si deve calcolare il doppio dell'area: AREA2

      AREA2=XJ*YK+XK*YI+XI*YJ-XK*YJ-XI*YK-XJ*YI
      IF (AREA2.LT.1.E-4) THEN
      PRINT*,' area negativa degenere '
      STOP
      END IF
      DO 10, I10=1,3
         DO 20, I20=1,6
            B(I10,I20)=B(I10,I20)/AREA2
   20    CONTINUE
   10 CONTINUE
      RETURN
      END