Informatica per Scienze

Geologiche LT

a.a.2017-2018

Introduzione all’utilizzo di metodologie

informatiche nella Geologia

Docente: Prof. Carla Braitenberg,

Dipartimento Matematica e Geoscienze, Via Weiss 1, Università di Trieste

E-mail: berg@units.it Tel. 040 5582258

Esempio da prima:

• x=0:1:40;

• y=10*exp(-abs(x));

• figure

• subplot(2,1,1)

• plot(x,y)

• subplot(2,1,2)

• semilogy(x,y,'k')

0 5 10 15 20 25 30 35 400

2

4

6

8

10

0 5 10 15 20 25 30 35 4010

-20

10-10

100

1010

Formulazione matematica della funzione

sinusoidale

• Le applicazioni delle funzioni di seno e coseno sono

ampie- ragione per la quale le analizziamo in dettaglio.

neoscillaziodell' ampiezzaA

onda.d'numerodelvettore

ydirezioneinondad'numero

xdirezioneinondad'numero

)sin(),(

k

k

k

ykxkAyxf

y

x

yx

Scriviamo la formula matematica di un’onda

sinusoidale in 3D. I fronti d’onda sono lineari ed hanno

un orientamento ben definito rispetto all’asse x. La

direzione del fronte d’onda viene espresso dal vettore

d’onda. Implementiamo poi tale equazione in uno

script matlab.

Definizione della direzione del numero d'onda

3.14162

k

)arctan(tanx

y

x

y

k

k

k

k

La direzione del vettore del numero d'onda è ortogonale ai

fronti d'onda e viene definita dalle due componenti del

vettore kx e ky. L'angolo del vettore con l'asse x e’ pari a:

X

Y

k

kx

ky

22 yx kkkk

La lunghezza del vettore

numero d’onda e’ pari al modulo

del vettore.

La lunghezza d’onda della

sinusoide e’ pari a:

Definizione della due componenti del vettore

numero d'onda

)sin(

)cos(

kk

kk

y

x

Da prima abbiamo definito la lunghezza del vettore

n’umero d’onda:

X

Y

k

kx

ky

22 yx kkkk

Le due componenti sono

allora date dalla relazione:

Matlab: onda sinusoidale lineare (sinusoide.m)

• function sinusoide(lam,A,alf)

• % grafico di fronte d'onda lineare

• % eliminazione variabili e librerie create in

precedenza

• if nargin==0

• sinusoide(25,10,30)

• else

• % creazione vettore x di valori compresi tra

0 e 70 uniformemente divisi su 100 valori.

Per y intervallo da 0 a 90.

• x = linspace(0,70,100);

• y= linspace(0,90,100);

• % creazione matrice a 2 Dimensioni X e Y.

• [X,Y]= meshgrid(x,y);

• % numero d'onda

• k=(2*pi)/lam;

• alf1=alf*pi/180;

• kvec=k*[cos(alf1) sin(alf1)];

• kx=kvec(1); ky=kvec(2);

• % calcolo dei valori asse z

• Z = A*(sin(kx*X+ky*Y));

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

% visualizzazione 3D (azimuth 37.5,

elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');

end

end

Valori tipo:

lam=25;

A=10

alf=30

Matlab: somma tra due funzioni

function sovrapposizione(lam,A,alf,a,b,c,d)

% grafico di fronte d'onda lineare

if nargin==0

sovrapposizione(15,3,45,0.08,0.02,0,0.003)

else

% creazione vettore x di valori compresi

% tra 0 e 70 uniformemente divisi su 100

% valori. Per y intervallo da 0 a 90.

x = linspace(0,70,100);

y= linspace(0,90,100);

% creazione matrice a 2 Dimensioni X e Y.

[X,Y]= meshgrid(x,y);

% numero d'onda

k=(2*pi)/lam;

alf1=alf*pi/180;

kvec=k*[cos(alf1) sin(alf1)];

kx=kvec(1); ky=kvec(2);

% calcolo dei valori asse z per fronte d'onda

Z1 = A*(sin(kx*X+ky*Y));

% calcolo dei valori asse z per polinomiale

Z2 = (a*X+b*Y+c*X.^3+d*Y.^2);

% sommo fronte d'onda + polinomiale

Z=Z1+Z2;

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

% visualizzazione 3D

% (azimuth 37.5, elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');

end

Esempio oscillazione

Onda sferica

• L'onda sferica rappresenta un'immagine istantanea

delle onde che si sviluppano sulla superficie di un lago

sollecitata dalla caduta di un masso.

• Descrizione matematica dell'onda sferica:

)sin()(sin 22 yxkArkAz

Se l'origine dell'onda e' in x0,y0:

))()(sin( 2

0

2

0 yyxxkAz

Matlab: onda sinusoidale lineare (sinusoide.m)

• function sinusoide_radiale(lam,A,x0,y0)

• % grafico di fronte d'onda lineare

• % eliminazione variabili e librerie create in

precedenza

• if nargin==0

• sinusoide_radiale(25,10,50,50)

• else

• % creazione vettore x di valori compresi tra

0 e 70 uniformemente divisi su 100 valori.

Per y intervallo da 0 a 90.

• x = linspace(0,100,100);

• y= linspace(0,100,100);

• % creazione matrice a 2 Dimensioni X e Y.

• [X,Y]= meshgrid(x,y);

• % numero d'onda

• k=(2*pi)/lam;

•

• % calcolo dei valori asse z

• Z = A*sin(k*sqrt((X-x0).^2+(Y-y0).^2));

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

% visualizzazione 3D (azimuth 37.5,

elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');

end

Valori tipo:

lam=25;

A=10

X0=50

Y0=50

Esempio onda sferica

Onda sferica la cui ampiezza

decade esponenzialmente

• Vogliamo anche descrivere il caso in cui

l'ampiezza dell'onda decade con il raggio

esponenzialmente.

))()(

(

0

2

0

2

0

2

20

2

20

),(

))()(sin(),(

sigy

yy

sigx

xx

eAyxA

yyxxkyxAz

Matlab: onda sinusoidale radiale che decade

function sinusoide_radexp(lam,A,x0,y0,sigx,sigy)

% grafico di onda radiale che decade

% eliminazione variabili e librerie create in

precedenza

if nargin==0

sinusoide_radexp(25,10,50,50,20,40)

else

x = linspace(0,100,100);

y= linspace(0,100,100);

% creazione matrice a 2 Dimensioni X e Y.

[X,Y]= meshgrid(x,y);

% numero d'onda

k=(2*pi)/lam;

% calcolo Ampiezza onda

Ar = A*exp(-((X-x0)/sigx).^2 -((Y-y0)/sigy).^2);

% calcolo dei valori asse z

Z = Ar.*sin(k*sqrt((X-x0).^2+(Y-y0).^2));

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

% visualizzazione 3D (azimuth 37.5,

elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');

end

Valori tipo:

lam=25;

A=10

x0=50

y0=50

sigx=20

sigy=40

Onda sferica con decadimento anisotropo

0

20

40

60

80

100

0

20

40

60

80

100

-505

XY

Z

Superficie polinomiale

• Di larga applicazione e' l costruzione di una superficie

costituita elementi polinomiali.

• Esempio: descrizione di un campo regionale di una

quantita' misurata, come emissione Radon, flusso

termico.

Esempio superficie polinomiale

• Z = 2*X+3*Y+X.^3*.1+Y.^2 *1.5

• Z = 0.2*X+0.3*Y+X.^3*0.01+Y.^2*0.05

• Esercizio:

– Costruire la superficie polinomiale. Variare i

parametri e descrivere il risultato.

– Costruire una superficie piana.

– Costruire una superficie di secondo ordine

function Polinomiale(a,b,c,d)

% grafico di onda radiale che decade

% eliminazione variabili e librerie create in

precedenza

if nargin==0

Polinomiale(0.2,0.3,0.01,0.05)

else

% creazione vettore x di valori compresi

tra 0 e 70 uniformemente divisi su 100

valori. Per y intervallo da 0 a 90.

x=linspace(-7,7,100);

y=linspace(-9,9,100);

% creazione matrice a 2 Dimensioni X e Y.

[X,Y]= meshgrid(x,y);

% calcolo dei valori asse z

Z = a*X+b*Y+c*X.^3+d*Y.^2;

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

shading flat

% visualizzazione 3D (azimuth 37.5,

elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');

end

end

Funzione polinomiale

-5

0

5

-5

0

5

-5

0

5

10

XY

Z

Graficare il gradiente di una superficie

• v = -2:0.2:2;

• [x,y] = meshgrid(v);

• z = y.*x.*exp(-x.^2 - y.^2);

• % da provare anche:

• % z = x.*exp(-x.^2 - y.^2);

• % z = y.*x.*exp(-x.^2 - y.^2);

• [px,py] = gradient(z,.2,.2);

• contour(v,v,z)

• hold on

• quiver(v,v,px,py)

• hold off

• figure

• surf(x,y,z);

• % visualizzazione 3D (azimuth 37.5, elevation 30)

• view(3);

Gradiente della topografia%Plot_DTM_isolinea

Zmin=180; Zmax=500;

step=2;

v=Zmin:step:Zmax;

sv=size(v);

load DTM_ortom_reg_regrid10m_WGS84_LD.mat;

figure('name',['esempio 2D isolinee Zmin: ' num2str(Zmin) ' Zmax: ' num2str(Zmax)])

%contour rappresenta in grafico le isolinee in 2D

[C, h]=contour(X1,Y1,Z1,v);

%axis([xmin xmax ymin ymax])

axis([13.74 13.76 45.715 45.725]);

figure

[px,py] = gradient(Z1,1,1);

contour(X1,Y1,Z1)

hold on

quiver(X1,Y1,px,py,5)

axis([13.74 13.76 45.715 45.725]);

13.74 13.742 13.744 13.746 13.748 13.75 13.752 13.754 13.756 13.758 13.7645.715

45.716

45.717

45.718

45.719

45.72

45.721

45.722

45.723

45.724

45.725

13.74 13.742 13.744 13.746 13.748 13.75 13.752 13.754 13.756 13.758 13.7645.715

45.716

45.717

45.718

45.719

45.72

45.721

45.722

45.723

45.724

45.725

Ulteriori applicazioni

• Creare un grafico di due onde sinusoidali di lunghezza

d’onda diversa sovrapposte. La direzione delle due sia la

medesima. Scegliere 1>> 2

• Creare un grafico di una superficie polinomiale di

secondo ordine che rappresenta l’andamento regionale,

al quale viene sovrapposta una anomalia a forma

gaussiana, che rappresenta un disturbo locale. La

larghezza della Gaussiana e’ determinata dal

coefficiente in denominatore dell’esponenziale.

• Simulare la presenza di tre mud-volcano; un mud-

volcano viene rappresentato tramite una funzione

gaussiana.

Esempio sovrapposizione

close all;

clear variables;

% parametri polinomiale

a = 1e-1;

b = 1e-1;

c = 1e-6;

d = 1e-6;

% parametri sinusoide

lam = 25;

A = 10;

x0 = 50;

y0 = 50;

sigx = 30;

sigy = 35;

% linspace x,y

x=linspace(-50,150,600);

y=linspace(-50,150,600);

% creazione matrice a 2 Dimensioni X e Y.

[X,Y]= meshgrid(x,y);

% POLINOMIALE

% calcolo dei valori asse z

Z_pol = a*X+b*Y+c*X.^3+d*Y.^2;

% SENO

% numero d'onda

k=(2*pi)/lam;

% calcolo Ampiezza onda

Ar_1 = A*exp(-((X-x0)/sigx).^2 - ...

((Y-y0)/sigy).^2);

% calcolo dei valori asse z

Z_sin_1 = Ar_1.*sin(k*sqrt((X-x0).^2+ ...

(Y-y0).^2));

% calcolo somma

Z = Z_pol + Z_sin_1;

% grafico tridimensionale dei risultati

figure

surf(X,Y,Z);

shading flat

% visualizz. 3D (azimuth 37.5, elevation 30)

view(3);

% manteniamo il fattore di scala

axis equal;

% etichette sugli assi

xlabel('X');

ylabel('Y');

zlabel('Z');