Le figure non vengono visualizzate, nel caso chiedetemi via mail che ve le invio
1 Primo
elaborato
1.1 Problema
assegnato
Il primo
problema prevede la realizzazione di una function MATLAB basata sull’algoritmo
del metodo delle corde modificate. Questo algoritmo, data una funzione f(x)
continua in [a0,b0] t.c. f(a0)*f(b0)<0, calcola le radici reali della
funzioni. Si riporta in allegato l’algoritmo sopracitato.
Successivamente
si applica tale function per il calcolo della radice della funzione
nell’intervallo
[0,5] e quella della funzione :
Nella parte finale
della relazione si analizzano i risultati ottenuti, confrontandoli con la
funzione “fzero”.
1.2 Functions MATLAB
1.2.1. Function “grafico”
La function “grafico” permette di individuare
l’intervallo nel quale la funzione intercetta l’asse delle x. L’algoritmo per
questa function è basato sul teorema degli zeri.
File:grafico.m
function [ ] = grafico(
a,b,f )
% Questa function crea il diagramma di una
funzione e indica se vale il teorema
% degli zeri nell’intervallo [a,b]
if
feval(f,a)*feval(f,b)>0;
disp (sprintf('\n
f(a)*f(b)>0 ->in questo intervallo non vale il teorema degli zeri'))
else
disp
(sprintf('\n
f(a)*f(b)<0 ->in questo intervallo vale il teorema degli zeri'))
end
% Finestra
p=a:0.01:b;
y=feval(f,p);
plot(p,y)
title ('Grafico funzione')
xlabel('x')
ylabel('y')
grid on
end
1.2.2. Function “cordemodificato”
La function
“cordemodificato” realizza l’algoritmo riguardante il metodo delle corde
modificate.
File:cordemodificato.m
function [ zo,fzo,n,ind
] = cordemodificato( a0,b0,f,tollf,Nmax)
% Metodo delle corde modificato basato sul
teorema degli zeri.
% L'obiettivo è calcolare la radice della
funzione f(x) con f:R->R
% INPUT:a0 limite sinistro dell'intervallo, b0
limite destro
% dell'intervallo, f funzione della quale si vuole
calcolare la radice,
% tollf tolleranza per il criterio di arresto
sulla funzione, Nmax è il
% numero max di iterazioni eseguibili.
% OUTPUT: zo zero della funzione, fzo valore
f(z), n numero di iterazioni,
% ind indicatore dell'errore: 0 criterio di
arresto verificato, 1 numero
% max di iterazioni raggiunto, -1 non
verificata la condizione di cambio di
% segno del teorema degli zeri.
F=feval(f,a0);
G=feval(f,b0);
if F*G>0;
ind=-1;
zo=[];
fzo=[];
n=[];
disp(sprintf('\n
f(a0)*f(b0)>0 non è valido il teorema degli zeri'))
return
end
w(1)=a0;
a(1)=a0;
b(1)=b0;
for n=1:Nmax
w(n+1)=(G*a(n)-F*b(n))/(G-F);
zo=w(n+1);
fzo=feval(f,w(n+1));
if feval(f,a(n))*fzo <= 0
a(n+1)= a(n);
b(n+1) = w(n+1);
G = fzo;
if feval(f,w(n))*fzo>0;
F = F/2;
end
else
a(n+1) = w (n+1);
b(n+1) = b(n);
F = fzo;
if feval (f,w(n))*fzo > 0
G = G/2;
end
end
if abs(fzo)<tollf;
fprintf('\n')
disp(' n zo f(zo)')
disp(sprintf(' %8.0f
%16.8e %16.8e', n, zo, fzo))
ind=0;
return
end
fprintf('\n')
disp(' n zo f(zo)')
disp(sprintf(' %8.0f %16.8e
%16.8e', n, zo, fzo))
ind=1;
end
1.2.3. Function “confronto”
Questa function
confronta i valori ottenuti dalla function “cordemodificato” con quelli della funzione
MATLAB “fzero”.
File:confronto.m
function [
err,errf,esito ] = confronto( f,c,zo,fzo )
%questa function confronta i valori ottenuti
con il metodo delle corde(function %cordemodificate) e i valori ricavati dalla
funzione fzero. INPUT: f funzione
%di cui si vuol calcolare la radice, c valore
vicino allo zero della
%funzione, zo valore dello zero della funzione
ottenuto con il metodo delle %corde, fzo valore della funzione in zo. OUTPUT:
err errore tra x e zo, errf
%errore tra f(x) e f(zo).
% Esito:
% esito=1 la function fzero ha trovato lo zero
della funzione, esito=-1
% l'algoritmo è stato terminato con output,
esito=-3 trovato valore NaN o
% Inf in un intervallo nel quale ho cambio di
segno, esito=-4 è stato
% trovato un valore nel campo complesso dove
la funzione cambia
% segno,esito=-5 fzero potrebbe convergere ad
un punto singolare, esito=-6
% fzero non rivela cambio di segno.
% Il formato [x,fval,exitflag]=fzero(f,c)
restituisce rispettivamente lo
% zero della funzione x, il valore della funzione
f(x) ed exitflag la
% condizione di uscita di fzero,
corrispondente in questo caso all'output
% esito.
[x,fval,exitflag]=fzero(f,c);
err=abs(x-zo);
errf=abs(fval-fzo);
esito=exitflag;
fprintf('\n')
disp(' x zo err f(x) f(zo) errf')
disp(sprintf('%16.8e %16.8e %16.8e
%16.8e %16.8e %16.8e',x,zo,err,fval,fzo,errf))
end
1.2.4. Function
“cub” e la function
“esp”
Le seguenti
function definiscono le funzioni che verranno utilizzate nell’elaborato.
File:cub.m
function [ y ] = cub( x
)
%Funzione x.^3-x-1 della quale dobbiamo
calcolare le radici reali
y = x.^3-x-1;
end
File:esp.m
function [ y ] = esp( x
)
%Funzione
y = exp(-4*x)-(1/10)
y = exp(-4*(x))-(1/10);
end
1.3 Applicazione 1
Calcolo della
radice reale di
nell’intervallo [0,5].
>>
grafico(0,5,@esp)
f(a)*f(b)<0 ->in questo intervallo vale
il teorema degli zeri
>> [zo fzo
n]=cordemodificato(0,5,@esp,eps,200)
n zo f(zo)
1
4.50000001e+000 -9.99999848e-002
n zo f(zo)
2
4.05000007e+000 -9.99999079e-002
n
zo f(zo)
3
3.31363698e+000 -9.99982476e-002
n zo f(zo)
4
2.29406874e+000 -9.98965348e-002
n zo f(zo)
5
1.21509860e+000 -9.22525724e-002
n zo f(zo)
6
4.60256650e-001 5.86544679e-002
n zo f(zo)
7
7.53648218e-001 -5.09341925e-002
n zo f(zo)
8
6.17286823e-001 -1.53429840e-002
n zo f(zo)
9
5.63351303e-001 5.04092829e-003
n zo f(zo)
10
5.76689522e-001 -4.16430014e-004
n zo f(zo)
11
5.75671734e-001 -1.01837083e-005
n zo f(zo)
12
5.75622155e-001 9.64794671e-006
n zo f(zo)
13
5.75646274e-001 -4.91261101e-010
n zo f(zo)
14
5.75646273e-001 -2.37032616e-014
n zo f(zo)
15
5.75646273e-001 2.37310172e-014
n zo f(zo)
16
5.75646273e-001 1.38777878e-017
zo =
0.5756
fzo =
1.3878e-017
n =
16
>> [err errf
esito]=confronto(@esp,0.5756,zo,fzo)
x zo err f(x) f(zo) errf
5.75646273e-001 5.75646273e-001 0.00000000e+000 1.38777878e-017 1.38777878e-017 0.00000000e+000
err =
0
errf =
0
esito =
1
1.4 Applicazione 2
Calcolo della
radice reale di
nell’intervallo [-3,3].
>>
grafico(-3,3,@cub)
f(a)*f(b)<0 ->in questo intervallo vale
il teorema degli zeri
>> [zo fzo
n]=cordemodificato(1,2,@cubica,eps,200)
n zo f(zo)
1
1.16666667e+000 -5.78703704e-001
n zo f(zo)
2
1.32330827e+000 -6.00390115e-003
n zo f(zo)
3
1.32654297e+000 7.79623624e-003
n zo f(zo)
4
1.32471556e+000 -1.02213517e-005
n zo f(zo)
5
1.32471795e+000 -1.73619750e-008
n zo f(zo)
6
1.32471796e+000 1.73029158e-008
n zo f(zo)
7
1.32471796e+000 2.22044605e-016
n zo f(zo)
8
1.32471796e+000 -8.88178420e-016
n zo f(zo)
9
1.32471796e+000 2.22044605e-016
n zo f(zo)
10
1.32471796e+000 2.22044605e-016
n zo f(zo)
11
1.32471796e+000 2.22044605e-016
n zo f(zo)
12
1.32471796e+000 -8.88178420e-016
n zo f(zo)
13
1.32471796e+000 2.22044605e-016
n zo f(zo)
14
1.32471796e+000 2.22044605e-016
n zo f(zo)
15
1.32471796e+000 2.22044605e-016
zo =
1.3247
fzo =
2.2204e-016
n =
15
>> [err errf
esito]=confronto(@cubica,1.32,zo,fzo)
x zo err f(x) f(zo) errf
1.32471796e+000 1.32471796e+000 0.00000000e+000 2.22044605e-016 2.22044605e-016 0.00000000e+000
err =
0
errf =
0
esito =
1
2 Secondo elaborato
Collasso per instabilità di una trave semplicemente appoggiata
2.1 Problema
Si
definisce carico critico,
per la teoria elastica della trave,
quella forza di compressione il cui valore porta indefinitamente ad inflessione
il solido snello su cui agisce, generando instabilità per carico di punta.
In questa relazione verrà trattato questo problema di stabilità delle
strutture, che tra i fenomeni di instabilità è quello più importante.
L’equazione
differenziale del 4° ordine della linea elastica
per una
generica trave inflessa (si definiscono inflesse le travi piane ad
asse rettilineo soggette a soli carichi flettenti, ossia normali all'asse della
trave) è:
Studiamo
il caso di trave prismatica (EJ=cost).
Valutate le quattro condizioni al contorno (date dalle condizioni statiche e da
quelle cinematiche), la risoluzione dell’equazione differenziale può essere
effettuata in termini di v(z).
Dopo aver determinato v(z), per derivazioni
successive, possono essere ricavati i valori della rotazione, della curvatura ,
del momento e del taglio attraverso le seguenti relazioni:
Lo studio delle travi si è sviluppato facendo l’ipotesi di
piccoli spostamenti:
(1) linearizzazione
delle relazioni cinematiche;
(2) Equilibrio
con riferimento alla configurazione in deformata.
Tale
ipotesi, senz’altro accettabile per le travi soggette a carichi trasversali,
non è più veritiera nel caso di carichi agenti lungo l’asse della trave,
soprattutto di compressione. In questo caso dobbiamo considerare anche la
presenza di uno sforzo assiale P.
Rimuovendo la seconda ipotesi, l’equazione della la linea elastica
assume questa forma:
Consideriamo il caso di assenza di carichi trasversali (q(z)=0).
L’equazione diventa:
Essendo
ottengo
;
Questa
è l’equazione che tratteremo relativa alla trave semplicemente appoggiata
soggetta a carico di punta P, schematizzata nella figura seguente:
Le
condizioni al contorno cinematiche sono
mentre le due condizioni statiche sono
Integrando l’equazione precedente ottengo:
Dividendo per
e ponendo
ottengo
La soluzione è del tipo:
con le seguenti condizioni al contorno:
(1)
(2)
Per la seconda condizione posso avere:
(1)
equilibrio
nella configurazione indeformata (configurazione rettilinea)
(2)
(altre
configurazioni di equilibrio)
Il
valore più piccolo di n (
) corrisponde al carico critico, pari quindi a
.La linea
elastica assume la seguente forma:
2.1 La modellazione del problema
Per quanto riguarda la discretizzazione
utilizzo il metodo di approssimazione alle differenze finite.
Dopo aver discretizzato il problema e aver
ottenuto un problema lineare del tipo
(dove
), determino il più piccolo autovalore che è quello relativo al carico critico
minimo P. Dopo aver definito
l’autovettore relativo a tale autovalore confronto i risultati ottenuti con il
procedimento analitico.
Per approssimare utilizzo la formula alle differenze
centrate:
Il problema lineare dato
,
,
assume la
forma
Si
può dimostrare che tale equazione corrisponde a
, dove
è
l’autovettore corrispondente a
e A la
matrice dei coefficienti che assume la seguente forma:
Nell’elaborato
successivo dovrò quindi ricercare gli autovalori e gli autovettori della matrice
A.
2.2 Function “car_critico”
Questa function ha come
obiettivo quello di determinare il carico critico di una trave.
File:car_critico.m
function [t,lam,Pc] =
car_critico(l,E,B,H,h)
% Questa Function
calcola il carico critico di una
trave
%appoggio-appoggio soggetta a una forza assiale P
%INPUTS
%Gli input riguardano le caratteristiche geometriche della
trave e l’intervallo
%di suddivisione di quest’ultima
% l[cm]=luce della trave
% E[KN/cm2]=modulo di elasticità del materiale
% B[cm]=larghezza della sezione
% H[cm]=altezza della sezione
% h=valore degli intervalli che suddividono la trave in n
parti
%OUTPUTS
% t=autovettore relativo al carico critico
% lam=autovalore minimo
% Pc[KN]=carico critico euleriano
n=l/h; % n° delle parti nelle quali è stata suddivisa la
trave
% Per la matrice A utilizzo
il metodo alle differenze finite centrate
A=diag(-2*ones(n,1))+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1);
% La funzione eig produce una matrice X avente per colonne
gli autovettori
%della matrice A, ed una matrice D diagonale avente sulla
stessa gli autovalori
%della matrice oggetto di studio;
% NB: all’ i-esimo autovalore corrisponde l'autovettore
associato
% alla colonna i-esima
[X,D]=eig(A);
lam=-(D(1,1));
s=1;
for i=2:n
if -(D(i,i))<lam;
lam=-(D(i,i));
s=i;
end
end
t=(X(:,s));
% autovettore associato a lam (autovalore
minimo)
J=B*H^3/12;
% momento di inerzia della sezione
Pc=lam*E*J/h.^2;
% carico critico di instabilità
End
2.2 Confronto dei risultati ottenuti
Esempi
per la validazione dell’algoritmo.
Dati dell’ esempio 1:
|
IMPUT
|
|
|
E (KN/cm2)
|
2.1
|
|
H( cm)
|
20
|
|
B (cm)
|
25
|
|
L (cm)
|
200
|
Al variare di h ottengo:
|
CAS0
|
h
|
Pc
|
|
1
|
20
|
7.08872 e+001
|
|
2
|
10
|
7.81842 e+001
|
|
3
|
5
|
8.21575 e+001
|
|
4
|
2
|
8.46505 e+001
|
|
5
|
1
|
8.55001 e+001
|
|
6
|
0.1
|
8.59284 e+001
|
SOLUZIONE ANALITICA
Dati dell’esempio 2:
|
IMPUT
|
|
|
E (KN/cm2)
|
2.1
|
|
H( cm)
|
30
|
|
B (cm)
|
40
|
|
L (cm)
|
150
|
Al variare di h ottengo:
|
CAS0
|
h
|
Pc
|
|
1
|
10
|
7.26316 e+001
|
|
2
|
5
|
7.75758 e+001
|
|
3
|
1
|
8.18072 e+001
|
|
4
|
0.5
|
8.23539 e+001
|
|
5
|
0.1
|
8.27942 e+001
|
SOLUZIONE ANALITICA
Possiamo
notare come i risultati ottenuti con il problema discreto convergano verso il
risultato conseguito per via analitica; infatti al diminuire del passo di
discretizzazione il risultato si avvicina sempre di più alla soluzione
analitica.
Nessun commento:
Posta un commento