Esercitazioni MatLab per l'esame di Metodi Numerici Ing. Navale a.a. 23-24

Y=Y(:)';

Z=Z(:)';

V=[X;Y;Z]; %matrice con 3 righe, coordinate X, Y e Z

subplot(1,2,1)

plottavettori(V,[0 0.4470 0.7410],0.1) %grafico vettori con applicazione in 0 e

componenti [tripletta RGB per colore]

axis equal

hold on

B=eye(3); % matrice identità

plottavettori(B,'r',2)

title('Vettori e base non trasformati');

legend('Vettori sfera','Base canonica',Location='southoutside');

xlabel('X');

ylabel('Y');

zlabel('Z');

rango=input('Inserire il rango della matrice A, da 1 a 3: '); % modifico rango per

avere ellissoide, ellisse o retta

if rango==1

A=[1 2 3;3 6 9;2 4 6];

elseif rango==2

A=[1 7 8;3 4 7;2 8 10];

else A=[1 7 0;3 4 0;2 8 2];

end

W=A*V; %moltiplico vettori per A

subplot(1,2,2)

plottavettori(W,[0 0.4470 0.7410],0.5)

axis equal

hold on

D=A*B; % moltiplico base canonica

plottavettori(D,'r',2)

hold on

[Autov,Autovett]=eig(A);

plottavettori(Autovett,'g',3) % plotto autovettori

title('Vettori e base trasformati');

legend('Vettori sfera','Base canonica','Autovettori',Location='southoutside');

xlabel('X');

ylabel('Y');

zlabel('Z');

str=sprintf('Rango di A=%d',rango);

text(0,3,14,str);

function plottavettori(V,colore,spessore)

zero=zeros(size(V));

X1=zero(1,:);

X=V(1,:);

Y=V(2,:);

Z=V(3,:);

quiver3(X1,X1,X1,X,Y,Z,'Color',colore,'LineWidth',spessore);

end

Codice per applicazione

classdef applicazione < matlab.apps.AppBase

% Properties that correspond to app components

properties (Access = public)

UIFigure matlab.ui.Figure

RangomatriceAKnob matlab.ui.control.DiscreteKnob

RangomatriceAKnobLabel matlab.ui.control.Label

FaccesferaunitariaEditField matlab.ui.control.NumericEditField

FaccesferaunitariaEditFieldLabel matlab.ui.control.Label

RunButton matlab.ui.control.Button

UIAxes2 matlab.ui.control.UIAxes

UIAxes matlab.ui.control.UIAxes

end

methods (Access = private)

function plottavettori_app(app,assi,V,colore,spessore) % funzione per plottare

zero=zeros(size(V));

X1=zero(1,:);

X=V(1,:);

Y=V(2,:);

Z=V(3,:);

quiver3(assi,X1,X1,X1,X,Y,Z,'Color',colore,'LineWidth',spessore);

axis (app.UIAxes, "equal")

end

end

% Callbacks that handle component events

methods (Access = private)

% Button pushed function: RunButton

function RunButtonPushed(app, event)

cla(app.UIAxes); % pulisce gli assi

cla(app.UIAxes2);

N=app.FaccesferaunitariaEditField.Value;

sphere(N); % grafico la sfera di NxN facce

[X,Y,Z]=sphere(N);

close all

X=X(:)'; % metto una sotto l'altra le colonne e traspongo, diventa una riga

Y=Y(:)';

Z=Z(:)';

V=[X;Y;Z]; % matrice di 3 righe, ogni riga sono le coordinate X, Y e Z

plottavettori_app(app,app.UIAxes,V,'b',0.8)

grid(app.UIAxes,"on")

hold(app.UIAxes,"on")

B=eye(3); % base canonica

plottavettori_app(app,app.UIAxes,B,'r',2)

legend(app.UIAxes,'Vettori sfera','Base canonica',Location='southeastoutside');

rango=app.RangomatriceAKnob.Value;

rango=str2double(rango); % da char faccio diventare numero (double) altrimenti non va

if

if rango==1

A=[1 2 3;3 6 9;2 4 6];

elseif rango==2

A=[1 7 8;3 4 7;2 8 10];

elseif rango==3

A=[1 7 0;3 4 0;2 8 2];

end

W=A*V;

plottavettori_app(app,app.UIAxes2,W,'b',0.8)

grid(app.UIAxes2,"on")

hold(app.UIAxes2,"on")

D=A*B;

plottavettori_app(app,app.UIAxes2,D,'r',2)

hold(app.UIAxes2,"on")

[Autov,Autovett]=eig(A);

plottavettori_app(app,app.UIAxes2,Autovett,'g',3)

legend(app.UIAxes2,'Vettori sfera','Base

canonica','Autovettori',Location='southeastoutside');

str=sprintf('Rango di A=%d',rango); % stampo il rango di A

text(app.UIAxes2,3,14,str);

end

% Value changed function: RangomatriceAKnob

function RangomatriceAKnobValueChanged(app, event)

changingValue = event.Value; %pulisce assi se inserisco nuovo valore del

knob cla(app.UIAxes);

cla(app.UIAxes2);

end

% Value changed function: FaccesferaunitariaEditField

function FaccesferaunitariaEditFieldValueChanged(app, event)

changingValue = event.Value; %pulisce assi se inserisco nuovo valore

dell'edit field

cla(app.UIAxes);

cla(app.UIAxes2);

end

end

% Component initialization

methods (Access = private)

% Create UIFigure and components

function createComponents(app)

% Create UIFigure and hide until all components are created

app.UIFigure = uifigure('Visible', 'off');

app.UIFigure.Position = [100 100 857 567];

app.UIFigure.Name = 'MATLAB App';

% Create UIAxes

app.UIAxes = uiaxes(app.UIFigure);

title(app.UIAxes, 'Vettori e base non trasformati')

xlabel(app.UIAxes, 'X')

ylabel(app.UIAxes, 'Y')

zlabel(app.UIAxes, 'Z')

app.UIAxes.Position = [343 298 414 259];

% Create UIAxes2

app.UIAxes2 = uiaxes(app.UIFigure);

title(app.UIAxes2, 'Vettori e base trasformati')

xlabel(app.UIAxes2, 'X')

ylabel(app.UIAxes2, 'Y')

zlabel(app.UIAxes2, 'Z')

app.UIAxes2.Position = [343 35 414 264];

% Create RunButton

app.RunButton = uibutton(app.UIFigure, 'push');

app.RunButton.ButtonPushedFcn = createCallbackFcn(app, @RunButtonPushed,

true); app.RunButton.Position = [125 190 100 23];

app.RunButton.Text = 'Run';

% Create FaccesferaunitariaEditFieldLabel

app.FaccesferaunitariaEditFieldLabel = uilabel(app.UIFigure);

app.FaccesferaunitariaEditFieldLabel.HorizontalAlignment = 'right';

app.FaccesferaunitariaEditFieldLabel.Position = [61 483 110 22];

app.FaccesferaunitariaEditFieldLabel.Text = 'Facce sfera unitaria';

% Create FaccesferaunitariaEditField

app.FaccesferaunitariaEditField = uieditfield(app.UIFigure, 'numeric');

app.FaccesferaunitariaEditField.ValueChangedFcn = createCallbackFcn(app,

@FaccesferaunitariaEditFieldValueChanged, true);

app.FaccesferaunitariaEditField.Position = [186 483 100 22];

% Create RangomatriceAKnobLabel

app.RangomatriceAKnobLabel = uilabel(app.UIFigure);

app.RangomatriceAKnobLabel.HorizontalAlignment = 'center';

app.RangomatriceAKnobLabel.Position = [133 280 94 22];

app.RangomatriceAKnobLabel.Text = 'Rango matrice A';

% Create RangomatriceAKnob

app.RangomatriceAKnob = uiknob(app.UIFigure, 'discrete');

app.RangomatriceAKnob.Items = {'1', '2', '3'};

app.RangomatriceAKnob.ValueChangedFcn = createCallbackFcn(app,

@RangomatriceAKnobValueChanged, true);

app.RangomatriceAKnob.Position = [143 317 60 60];

app.RangomatriceAKnob.Value = '3';

% Show the figure after all components are created

app.UIFigure.Visible = 'on';

end

end

% App creation and deletion

methods (Access = public)

% Construct app

function app = applicazione

% Create UIFigure and components

createComponents(app)

% Register the app with App Designer

registerApp(app, app.UIFigure)

if nargout == 0

clear app

end

end

% Code that executes before app deletion

function delete(app)

% Delete UIFigure when app is deleted

delete(app.UIFigure)

end

end

end Esercitazione 3: sistema massa-molla-smorzatore

Vogliamo simulare il sistema massa-molla-smorzatore (vedere le dispense di Bruzzone):

( F

c k

ẍ(t) + ẋ(t) + x(t) = sin(ωt)

0

m m m

x(0) = x , ẋ(0) = v

0 0

p

dove c = 2 ζ ω m e ω = k/m.

n n

Prima di applicare un metodo numerico, occorre riscrivere l’equazione differenziale del secondo

ordine come un sistema di due equazioni differenziali del primo ordine:

 ẋ(t) = v



 F

k

c v x + sin(ωt)

v̇(t) = 0

−

− m m m

 x(0) = x , v(0) = v

 0 0

Implementare il metodo Runge-Kutta di ordine 4 (funzione di Matlab ode45) e il metodo di

Eulero Esplicito (ordine 1).

Usare i seguenti valori: posizione e velocità iniziali nulle, tempo finale 15 secondi, m = 20,

k = 20, F = 5, ω = 2, ζ = 2 oppure ζ = 0.2.

0

N.B Se volete vedere lo smorzamento, inserite F 0 = 0, altrimenti con una forzante sinusoidale

non si smorza mai.

clc

clear all

close all

metodo=input('Inserire metodo; 1 per Runge-Kutta (Ode45), 2 per Eulero esplicito: ');

% chiede metodo da usare

m=20; % massa

k=20; % costante elatsica

F0=5; % forzante

omega=2; % pulsazione

zita=2; % fattore di smorzamento (>1 sovrasmorzato, <1 sottosmorzato)

t0=0; % tempo iniziale

tf=15; % tempo finale

h=0.1; % passo temporale

y0=[0;0]; % condizioni iniziali

if metodo==1 % metodo Runge_Kutta (ode45)

[t,y]=ode45(@(t,y)fun(t,y,m,k,F0,omega,zita),[t0 tf],y0); % chiama la funzione

figure(1)

plot(t,y(:,1),'r-',t,y(:,2),'b-',Marker='o')

title('Soluzione massa-molla-smorzatore con CI nulle ODE45')

legend('Posizione x(t)','Velocità v(t)',Location='southeast')

xlabel('Tempo [s]')

ylabel('x(t),v(t)')

grid on

else % metodo Eulero esplicito

N=(tf-t0)/h; % numero passi temporali

y=zeros(N+1,2); % inizializzo matrice (N+1)x2 delle soluzioni

t=linspace(t0,tf,N+1); % vettore tempi

y(1,:)=y0'; % condizioni iniziali nella matrice soluzione

for n=1:N % ciclo for per calcorae soluzioni usando la precedente

y(n+1,:)=y(n,:)+h*eulero_esplicito(t(n),y(n,:),m,k,F0,omega,zita)'; % chiamo la

fun

funzione f

end figure(2)

plot(t,y(:,1),'r-',t,y(:,2),'b-',Marker='o')

title('Soluzione massa-molla-smorzatore con CI nulle Eulero esplicito')

legend('Posizioe','Velocità',Location='southeast')

xlabel('Tempo [s]')

ylabel('x(t),v(t)')

grid on

end

function f=fun(t,y,m,k,F0,omega,zita)

omega_n=sqrt(k/m); % pulsazione naturale sistema

c=2*zita*omega_n*m;

f=[y(2);F0/m*sin(omega*t)-c/m*y(2)-k/m*y(1)]; % vettore colonna che contiene derivate

y'(1) e y'(2)

end

function g=eulero_esplicito(t,y,m,k,F0,omega,zita) % funzione f

omega_n=sqrt(k/m);

c=2*zita*omega_n*m;

g=[y(2); F0/m*sin(omega*t)-c/m*y(2)-k/m*y(1)];

end Esercitazione Fast Fourier Transform

Lo scopo dell’esercitazione è effettuare la FFT di un segnale contenuto nel file ”signal1.mat” e dello

stesso segnale con noise contenuto nel file &