Matlab , newton raphson , quadrilatero articolato

di
jorian69
il
4 risposte
J
jorian69
Pochi posts
Iscritto da
dic, 2011
Messaggi:
5

Matlab , newton raphson , quadrilatero articolato

Buon Giorno e quasi Buon anno a tutti!
premetto che sono un profano nella progettazione di matlab.
ho riscontrato dei problemi nel trovare i parametri di un quadrilatero articolato. ho creato il codice, e non mi da problemi nel calcolo degli angoli della biella e della manovella. il problema inizia a crearsi quando provo a calcolare le velocità angolari, cosa che non riesco a spiegarmi.
NB : l' <x0_p> introdotto dovrebbe almeno restituire i valori di omega3 e omega4 dell'angolo 0°
visto che quel valore me lo sono calcolato con il metodo delle polari ( con matita e righello ).

vi allego la funzione principale e a seguire tutte le altre richiamate :
una tabella delle posizioni fi3 e fi4 per una sola configurazione possibile,
una tabella velocità omega3 e omega4, da dove però non riesce a stampare niente)
un grafico delle posizioni degli angoli fi3 e fi4 in funzione di fi2 (manovella)

function tabulaB( )
%calcolo delle posizioni di un quadrilatero articolato
caso =2;

x0 = [ .5 4.5 ] ;
x01 = [ 1 5.5 ] ;
tn = (pi/8)*(0:16);

col = zeros (length(tn),2);

for i=1:length(tn)
sol = NewRap( x0 , tn(i), caso );
solgrad = sol*180/pi ;
col(i,:) = RendiLegg ( solgrad , 360 ) ;
%--------------------------------------- primo controllo
if (tn(i)>pi)
sol = NewRap( x01 , tn(i), caso );
solgrad = sol*180/pi ;
col(i,:) = RendiLegg ( solgrad , 360 ) ;
end
%--------------------------------------fine primo controllo
%--------------ora creo i valori della configurazione alternativa

colla = [col(i,:),det(jacob(sol,caso))] ;
fprintf ('vars = %3.1f -----> ',tn(i)*180/pi) ; printa(colla);
end
%-------------------------------------tabella 2 : velocità
col_rad = col*pi/180;
col_p = zeros(length(tn),2);
tn_p = [ tn' omega2(tn')];
x0_p = [ 17.5 18.75 ];

for i=1:length(tn)
sol_p = NewRap_p ( x0_p , tn_p , col_rad(i,:) , caso );
sol_p = RendiLegg (sol_p,2*pi);
col_p(i,:) = sol_p;
fprintf ('vars = %3.1f ---%6.3f--> ',tn(i)*180/pi,omega2(tn(i))) ; printa(col_p(i,:));

end

%--------------------------------grafico------------------------
X = tn;
Y = col;
plot(X,Y);
legend ('fi3','fi4','location','NorthEastOutSide')

end



function [ var ] = funz( x , tn , caso )
%legge in ingresso il vettore delle <x> e di <tn> (termini noti)
% restituisce il vettore delle <var>,
%della stessa lunghezza del vettore delle <x>.
%sintassi :
% [var] = funz ( x , tn , caso)
%
% caso : swicha tra le tre funzioni:
% caso 1 -- manovellismo di spinta
% caso 2 -- quadrilatero articolato
% caso 3 -- guida di fairbairn

var = zeros (1,length(x));

switch caso

case 1
% esercitazione 4 A
var(1) = 0.2*cos(tn(1))+0.6*cos(x(1))-x(2) ;
var(2) =0.2*sin(tn(1))+0.6*sin(x(1)) ;
case 2
% esercitazione 4 B
var(1) = -0.7+0.3*cos(tn)+0.4*cos(x(1))+0.6*cos(x(2)) ;
var(2) = +0.3*sin(tn)+0.4*sin(x(1))+0.6*sin(x(2)) ;
case 3
%esercitazione 4 C - guida di Fairbairn
var(1) = 11.45*cos(tn)+x(1)*cos(x(2)+pi) ;
var(2) = 11.45*sin(tn)+x(1)*sin(x(2)+pi)+45.8 ;
var(3) = 80*cos(x(2))+x(4) ;
var(4) = 80*sin(x(2))+x(3)-100 ;
otherwise
end
end

function [ soluzioni ] = funzPrim ( x , tn , valFunz, caso )
% ESERCITAZIONE 5
%sintazssi : [sol] = funzPrim ( x , tn , valFunz , caso)
%
% derivate prime di funz
%
% x = vettore delle incognite
% tn = termine noto : [ fi2 omega2 ]
% omega2 = velocità angolare della manovella 2
% valFunz = valori trovati nella funz ----> da inizializzare
% soluzioni = valori restituiti

%caso = 1 : manovellismo di spinta
% 2 : quadrilatero articolato
% 3 : guida di fairbairn

switch caso
case 1
soluzioni = zeros(1,2);
% caso 1 : derivate prime di funz
% x = vettore delle incognite : [ omega3 esse_punto ]
omega3 = x(1);
esse_punto = x(2);
% tn = termini noti : [ fi2 omega2 ]
% omega2 = velocità angolare della manovella 2
fi2 = tn(1);
omega2 = tn(2);
% valFunz = valori trovati nella funz : [ fi3 esse theta]
fi3 = valFunz(1);

soluzioni(1) = -esse_punto -.2*omega2*sin(fi2) -.6*omega3*sin(fi3);
soluzioni(2) = 0 +.2*omega2*cos(fi2) +.6*omega3*cos(fi3);

case 2
soluzioni = zeros(1,2);
% caso 2 : derivate prime di funz
% x = vettore delle incognite : [ omega3 omega4 ]
omega3 = x(1);
omega4 = x(2);
% tn = termine noto : [ fi2 omega2 ]
% omega2 = velocità angolare della manovella 2
fi2 = tn(1);
omega2 = tn(2);
% valFunz = valori trovati nella funz : [ fi3 fi4 det(jacobiana) ]
fi3 = valFunz(1);
fi4 = valFunz(2);

soluzioni(1) = -.3*omega2*sin(fi2) -.4*omega3*sin(fi3) -.6*omega4*sin(fi4);
soluzioni(2) = +.3*omega2*cos(fi2) +.4*omega3*cos(fi3) +.6*omega4*cos(fi4);

case 3
soluzioni = zeros(1,4);
% caso 3 : derivate prime di funz
l2 = 11.45 ;
l4 = 80 ;
% x = vettore delle incognite : [ esse_p omega4 yB_p xC_p ]
esse_p = x(1);
omega4 = x(2);
yB_p = x(3);
xC_p = x(4);
% tn = termine noto : [ fi2 omega2 ]
% omega2 = velocità angolare della manovella 2
fi2 = tn(1);
omega2 = tn(2);
% valFunz = valori trovati nella funz : [ esse fi4 yB xC ]
esse = valFunz(1);
fi3 = valFunz(2);
yB = valFunz(3);
xC = valFunz(4);

soluzioni(1) = -l2*omega2*sin(fi2) +esse_p*cos(fi4+pi) -esse*omega4*sin(fi4+pi) ;
soluzioni(2) = +l2*omega2*sin(fi2) +esse_p*sin(fi4+pi) +esse*omega4*cos(fi4+pi) ;
soluzioni(3) = -l4*omega4*sin(fi4) -xC_p ;
soluzioni(4) = +l4*omega4*cos(fi4) +yB_p ;

otherwise
end

function [ Jacob ] = jacob( x , caso )
% legge il vettore incognite <x>
% e crea la matrice jacobiana <Jacob>, di dimensione (length(x),length(x))
% sintassi :
% [Jacob] = jacob ( x , caso )
%
%caso : swicha tra le tre funzioni:
% caso 1 -- manovellismo di spinta
% caso 2 -- quadrilatero articolato
% caso 3 -- guida di fairbairn

Jacob = zeros (length(x),length(x));

switch caso

case 1
% esercizio A
Jacob(1,1) = -0.6*sin(x(1)) ;
Jacob(1,2) = -1 ;
Jacob(2,1) = 0.6*cos(x(1));
Jacob(2,2) = 0 ;

case 2
% esercizio B
Jacob(1,1) = -0.4*sin(x(1)) ;
Jacob(1,2) = -0.6*sin(x(2)) ;
Jacob(2,1) = 0.4*cos(x(1));
Jacob(2,2) = 0.6*cos(x(2)) ;

case 3
% esercizio C - guida di fainberg
% x(1) = spostamento carrello glifo
% x(2) = angolo glifo
% x(3) = posizione su y
% x(4) = posizione su x
Jacob(1,1) = cos(x(2)+pi) ;
Jacob(1,2) = -x(1)*sin(x(2)+pi) ;
Jacob(2,1) = sin(x(2)+pi) ;
Jacob(2,2) = x(1)*cos(x(2)+pi) ;
Jacob(3,2) = -80*sin(x(2)) ;
Jacob(3,4) = 1 ;
Jacob(4,2) = 80*cos(x(2)) ;
Jacob(4,3) = 1 ;
otherwise
end

end

function [ X ] = itera( X0 , Jacob , Funz )
% passo elementare di iterazione
X = X0 - Funz*((inv(Jacob))') ;
end

function [ modul ] = modulo( x )
%ritorna il vettore <modulo>, composto dai moduli dei valori del vettore
%inserito <x>
N = length(x);
modul = zeros(1,N);
for i=1:N
if (x(i)>0)
modul(i) = x(i) ;
else
modul(i) = -x(i) ;
end
end
end

function [ Sol ] = NewRap( x0 , tn, caso )
%calcola il valore degli zeri della funzione studiata
% x0 : vettore valori di inizializzazione delle incognite
% tn : vettore valori termini noti
% caso : swicha tra le tre funzioni:
% caso 1 -- manovellismo di spinta
% caso 2 -- quadrilatero articolato
% caso 3 -- guida di fairbairn
% sintassi :

N = 200;
C = length(x0);
X = zeros (N,C);
X(1,:) = x0 ;
E = ones(1,C)*0.00001;
D = ones(1,C)*0.00001;
Sol = zeros(1,C);

% fprintf ('X0 =') ; printa (X(1,:)) ;

for i = 2:N

% fprintf (' F(%i) = ',i-2) ; printa (funz ( X(i-1,:), tn )) ;
% fprintf (' Jacobiano(%i) =',i-2) ; jacob( X(i-1,:) ) ;
% fprintf (' detJacob(%i) = %f \n',i-2,det(jacob( X(i-1,:) ))) ;
% fprintf ('\n InversaJacob(%i) = ',i-2) ; inv(jacob( X(i-1,:) ))
if (modulo(det(jacob( X(i-1,:),caso)))<1*10^(-6))
fprintf('errore nel determinante');
else
X(i,:) = itera ( X(i-1,:) , jacob(X(i-1,:),caso) , funz(X(i-1,:),tn,caso) );
% fprintf ('X(%i) = ', i-1 ) ; printa (X(i,:));
var1 = modulo(funz(X(i,:),tn,caso))<E ;
var2 = modulo( X(i,:)-X(i-1,:) )<D ;
if ( isequal( var1 & var2 , ones(1,C)) )
Sol = X(i,:);
break
end
end
end
end

function [ Sol ] = NewRap_p( x0 , tn , valFunz , caso )
%calcola il valore degli zeri della funzione studiata
% x0 : vettore valori di inizializzazione delle incognite
% tn : vettore valori termini noti
% valFunz : vettore valori trovati nella funz
% caso : swicha tra le tre funzioni:
% caso 1 -- manovellismo di spinta
% caso 2 -- quadrilatero articolato
% caso 3 -- guida di fairbairn
% sintassi : [ Sol ] = NewRap_p ( x0 , tn , valFunz , caso )


N = 2000;
C = length(x0);
X = zeros (N,C);
X(1,:) = x0 ;
E = ones(1,C)*0.001;
D = ones(1,C)*0.001;
Sol = zeros(1,C);


for i = 2:N

if (modulo(det(jacob( X(i-1,:),caso)))<1*10^(-6))
fprintf('errore nel determinante');
else
X(i,:) = itera ( X(i-1,:) , jacob(X(i-1,:),caso) , funzPrim(X(i-1,:),tn,valFunz,caso) );
var1 = modulo(funzPrim(X(i,:),tn,valFunz,caso))<E ;
var2 = modulo( X(i,:)-X(i-1,:) )<D ;
if ( isequal( var1 & var2 , ones(1,C)) )
Sol = X(i,:);
break
end

end
end
end

function [ Omega2 ] = omega2 ( fi2 )
% questa funzione calcola il valore di omega2 dell'ESERCITAZIONE 5
% SINTASSI: [omega2] = omega2 ( fi2 )

k = 0.012 ;
omega20 = 25 ; % rad/s

Omega2 = omega20*(1+k*sin(fi2)) ;
end

function [ ] = printa( x )
% stampa i valori del vettore x e va a capo

fprintf ('[')
for i = 1 : length(x)
fprintf('\t%f',x(i)) ;
end
fprintf ('\t]\n');
end

function [ a ] = RendiLegg( r , x )
% Rendi Leggibile.
% risolve l'equazione:
% r = a + k*x
%sintassi:
% [a] = RendiLegg ( r , x )

if (length(r)==1)
if ((r<=x/2) && (r>=-x/2))
% if ((r<=x) && (r>=0))
a = r ;
else if (r>0)
a = RendiLegg( r - x , x ) ;
else
a = RendiLegg( r + x , x ) ;
end
end
else
l = length(r) ;
a = zeros (1,l);
for i = 1:l
a(i) = RendiLegg ( r(i) , x ) ;
end

end

end

4 Risposte

  • J
    jorian69
    Pochi posts
    Iscritto da
    dic, 2011
    Messaggi:
    5

    Re: Matlab , newton raphson , quadrilatero articolato

    Se il codice vi risulta incomprensibile, posso aggiungere e migliorare i commenti. fatemi sapere per favore. grazie a tutti
  • G
    giug
    Molto attivo
    Iscritto da
    nov, 2011
    Messaggi:
    641

    Re: Matlab , newton raphson , quadrilatero articolato

    Purtroppo non è detto che l'argomento sia conosciuto da tutti, quindi avendo tutto quel codice davanti è anche difficile riuscire a capire quale sia la domanda.
    Io ho fatto girare il codice (la funzione principale) e arriva fino alla fine senza dare errori e facendo vedere il grafico. Non conosco l'argomento, quindi non so giudicare se i risultati siano coerenti o meno.
    Forse potresti provare a fare domande più precise, tipo "fino alla riga 10 non ci sono problemi, alla riga 11 dovrebbe dare questo risultato invece dà quell'altro".
  • C
    ciacco
    Pochi posts
    Iscritto da
    giu, 2015
    Messaggi:
    5

    Re: Matlab , newton raphson , quadrilatero articolato

    Mi interesserebbe capire il metodo Newton raphson come implementarlo e la logica che segue se fosse possibile al variare dell'angolo "ALFA" che si fornisce in ingresso per consentire il moto al quadrilatero, con seguente plottaggio delle posizioni dei link del quadrilatero e plot della traiettoria del punto C3 del quadrilatero in figura con posizione del telaio note

    Attualmente prendendo in considerazione questa parte di codice 
    function tabulaB( )
%calcolo delle posizioni di un quadrilatero articolato
caso =2;

x0 = [ .5 4.5 ] ;
x01 = [ 1 5.5 ] ;
tn = (pi/8)*(0:16);

col = zeros (length(tn),2);

for i=1:length(tn)
sol = NewRap( x0 , tn(i), caso );
solgrad = sol*180/pi ;
col(i,:) = RendiLegg ( solgrad , 360 ) ;
%--------------------------------------- primo controllo
if (tn(i)>pi)
sol = NewRap( x01 , tn(i), caso );
solgrad = sol*180/pi ;
col(i,:) = RendiLegg ( solgrad , 360 ) ;
end
%--------------------------------------fine primo controllo
%--------------ora creo i valori della configurazione alternativa

colla = [col(i,:),det(jacob(sol,caso))] ;
fprintf ('vars = %3.1f -----> ',tn(i)*180/pi) ; printa(colla);
end
%-------------------------------------tabella 2 : velocità
col_rad = col*pi/180;
col_p = zeros(length(tn),2);
tn_p = [ tn' omega2(tn')];
x0_p = [ 17.5 18.75 ];

for i=1:length(tn)
sol_p = NewRap_p ( x0_p , tn_p , col_rad(i,:) , caso );
sol_p = RendiLegg (sol_p,2*pi);
col_p(i,:) = sol_p;
fprintf ('vars = %3.1f ---%6.3f--> ',tn(i)*180/pi,omega2(tn(i))) ; printa(col_p(i,:));

end

%--------------------------------grafico------------------------
X = tn;
Y = col;
plot(X,Y);
legend ('fi3','fi4','location','NorthEastOutSide')

end
    mi da problemi sulla funzione NEWRAP
  • J
    jorian69
    Pochi posts
    Iscritto da
    dic, 2011
    Messaggi:
    5

    Re: Matlab , newton raphson , quadrilatero articolato

    Nel codice in alto ci sta solo il programma, mentre in basso ho richiamato tutte le funzioni utilizzate e create.

    Ti da problemi su NEWRAP perché ancora non esiste sul tuo MATLAB

    Ora, se ho capito bene la tua domanda, lo scopo di questo programma è proprio quello di trovare la posizione angolare del cedente in funzione della posizione angolare del movente, tuttavia la posizione trovata (esistono 2 zeri, e il metodo di newton-rapson si ferma quando ne trova uno) può non essere compatibile con il tipo di configurazione scelta (simmetrica o antisimmetrica, mi sembra di ricordare), quindi per verificare la veridicità dei risultati ti conviene usare un metodo grafico.

    Un'altra strada per trovare la soluzione è attraverso il metodo analitico: esiste la soluzione analitica, che puoi trovare in letteratura in molteplici libri, ma ora ricordo solo vagamente che stava su uno di quelli scritti anche da "Pennestrì" docente di Torvergata in "Meccanica Applicata alle Macchine". Ti consiglio la lettura di quasi tutta la sua collana, reperibile in molte biblioteche tecnologiche, mentre per aiuto diretto ora sono un po arruginito. Se stai a Roma possiamo incrociarci in un università e te li faccio vedere direttamente.
Devi accedere o registrarti per scrivere nel forum
4 risposte