Péndulo amortiguado

Un péndulo amortiguado sigue la ecuación:

$$ m \cdot l\cdot \theta''+k \cdot l\cdot \theta'+m \cdot g \cdot \sin(\theta)=0 $$

donde

  • $l$ longitud del péndulo (= 10 m).
  • $\theta$ ángulo que forma la cuerda con la vertical.
  • g aceleración de la gravedad (≈ 10 m/sg ).
  • m masa de la bola (= 20 kg).
  • k coeficiente de fricción del medio (= 4 kg/sg).

El péndulo comienza en la posición de reposo. Se imprime al péndulo una velocidad inicial de 1 rad/s. Calculamos la trayectoria del péndulo mediante del método de Euler.

Apartado 1

  • Resuelve de forma aproximada el PVI correspondiente al problema anterior, cuando y(0) = 1, usando un método adaptativo, en el intervalo de tiempo [0,60]. Dibuja los puntos que ha utilizado el método adaptativo para aproximar la solución.
  • Reescribimos la ecuación de segundo orden como un sistema de primer orden.
  • Usamos el método adaptativo Runge-Kutta 4,5, llamando a ode45. El argumento tspan es solo un array con el valor inicial y el final, en vez de un linspace, porque no queremos que ode45 haga la interpolación, sino que nos devuelva los puntos que ha usado para aproximar la solución de la ecuación diferencial.
In [1]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2); -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 30;
tspan = [t0, tf];
y0 = [0,1];

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

hold on;
plot(ts, thetas, 'b-')
plot(ts, thetas,'bo')
hold off;
xlabel('t')
ylabel('\theta')
title('Las oscilaciones se amortiguan...')

Apartado 2

  • Interpola los valores obtenidos en el método anterior usando el polinomio de Lagrange, y la spline cúbica. Estima el máximo de la diferencia entre el resultado de interpolar con uno y otro método, en el intervalo [0,60].
  • Usamos CubicSpline para la interpolación con splines.
In [2]:
% Evaluamos en un array de puntos x_eval, para dibujar el spline
t_eval = linspace(t0,tf,200);
y_eval = spline(ts, thetas, t_eval);
hold on;
plot(ts, thetas, 'o');
plot(t_eval, y_eval,'-');
title('Spline interpolador')
xlabel('x')
ylabel('y')
hold off;
legend('nodos de interpolacion', 'spline')
  • Para la interpolación con el polinomio de Lagrange, mostramos el resultado usando código explícito que usamos en uno de los cuadernos, y después polyfit, que es más estable.
In [3]:
function y_eval = lagrangepol(xs, ys, x_eval)
%    Polinomio k-esimo de la base de Lagrange
%    
%    INPUTS:
%        - `xs` coordenada x de los puntos de interpolación
%        - `ys` coordenada y de los puntos de interpolación
%        - `x_eval` es una coordenada, o un array de coordenadas x
%          donde evaluar el polinomio
%    
%    OUTPUTS
%        - `y_eval` array resultado de evaluar el polinomio interpolador
%        en x_eval
    n = length(xs);
    nx = length(x_eval);
    y_eval = zeros(1,nx);
    for i=1:n
        xi = xs(i);
        yi = ys(i);
        base = ones(1,nx);
        for j=1:n
            if j~=i
                xj = xs(j);
                base = base.*(x_eval - xj)./(xi - xj);
            end
        end
        y_eval = y_eval + base.*yi;
    end
end
In [4]:
% Evaluamos en un array de puntos x_eval, para dibujar el spline
t_eval = linspace(t0,tf,200);
y_eval = lagrangepol(ts, thetas, t_eval);
hold on;
plot(ts, thetas, 'o');
plot(t_eval, y_eval,'-');
title('Polinomio interpolador')
xlabel('x')
ylabel('y')
hold off;
legend('nodos de interpolacion', 'polinomio interpolador')
In [5]:
% Evaluamos en un array de puntos x_eval, para dibujar el polinomio interpolador
t_eval = linspace(t0,tf,200);
cs = polyfit(ts, thetas, length(ts)-1);
y_eval = polyval(cs, t_eval);
hold on;
plot(ts, thetas, 'o');
plot(t_eval, y_eval,'-');
title('Polinomio interpolador')
xlabel('x')
ylabel('y')
hold off;
legend('nodos de interpolacion', 'polinomio interpolador')

warning: matrix singular to machine precision, rcond = 1.78761e-70
warning: called from
    polyfit at line 119 column 5
  • Curiosidad: mostramos el resultado de pedir a ode45 que interpole: usa un polinomio de grado 4 a trozos, y es mucho mejor opción que cualquiera de las anteriores...
In [6]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2); -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 30;
y0 = [0,1];

N = 200;
tspan = linspace(t0, tf,N);

[ts_ode45, ys_ode45] = ode45(fun, [t0, tf], y0);

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

hold on;
plot(ts_ode45, ys_ode45(:,1), 'ro')
plot(ts, thetas, 'b-')
hold off;
xlabel('t')
ylabel('theta')
title('Las oscilaciones se amortiguan...')

Apartado 3

Fijamos los valores de

  • $l$ longitud del péndulo (= 10 m).
  • k coeficiente de fricción del medio (= 4 kg/sg).

pero buscamos elegir el parámetro m de modo que al dejar caer el péndulo desde la posición inicial $\theta(t_0)=-\pi/2$, el péndulo rebote hasta una posición máxima de $\theta(t_0)=+\pi/4$. Para ello:

  1. Escriba una función f que reciba como argumento la masa m del péndulo, y devuelva el ángulo máximo que alcanza el péndulo cuando parte de su posición inicial en reposo con $\theta(t_0)=-\pi/2$.
  2. Encuentre el valor de m para el que f(m) es igual a $\theta(t_0)=+\pi/4$. Explique qué método ha utilizado para encontrar ese valor de m.
  3. Compruebe la solución dibujando la trayectoria con las condiciones iniciales del enunciado, y la masa m encontrada.

1. Escriba una función posicion_final que reciba como argumento la masa m del péndulo, y devuelva el ángulo máximo que alcanza el péndulo cuando parte de su posición inicial en reposo con $\theta(t_0)=-\pi/2$.

In [7]:
function target = posicion_final(m)
    l = 10;
    k = 4;
    g = 10;
    %Usamos un intervalo temporal suficiente como para que el pendulo
    % haya tenido tiempo de rebotar
    t0 = 0;
    tf = 20;
    t_span = [t0,tf];
    y0 = [-pi/2,0];

    fun = @(t,y)([y(2); -(k/m)*y(2)-(g/l)*sin(y(1))]);

    beta = pi/4;
    [t,y] = ode45(fun, t_span, y0);
    ymax = max(y(:,1));
    target = ymax - beta;
end
In [8]:
posicion_final(1), posicion_final(20)

ans = -0.79605
ans =  0.24172
  • 2. Encuentre el valor de m para el que f(m) es igual a $\theta(t_0)=+\pi/4$. Explique qué método ha utilizado para encontrar ese valor de m.
  • Hemos usado el método de bisección, después de haber encontrado dos valores de la masa m que dan lugar a valores del rebote máximo que son, respectivamente, menor y mayor que el buscado.
In [9]:
m0 = fsolve(@posicion_final, 1)

m0 =  11.433

  • 3. Compruebe la solución dibujando la trayectoria con las condiciones iniciales del enunciado, y la masa m encontrada.

  • Dibujamos una línea horizontal a la altura que corresponde al ángulo buscado, para comprobar que el máximo rebote es exactamente el buscado.

In [10]:
m = m0;
l = 10;
k = 4;
g = 10;
beta = pi/4;

fun = @(t,y)([y(2); -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 30;
y0 = [-pi/2,0];

N = 200;
tspan = linspace(t0, tf,N);

[ts_ode45, ys_ode45] = ode45(fun, [t0, tf], y0);

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

hold on;
plot(ts_ode45, ys_ode45(:,1), 'ro')
plot(ts, thetas, 'b-')
plot([t0, tf], [beta, beta], 'g-')
hold off;
xlabel('t')
ylabel('\theta')
title('Las oscilaciones se amortiguan...')

Apartado 4

Volvemos a fijar la masa m=20.

Recordad que este modelo asume que la masa se concentra en el extremo del péndulo:

  • la velocidad de la masa situada en el extremo del péndulo es $l\cdot |\theta'(t)|$.
  • el desplazamiento $\Delta s$ es la distancia recorrida $\Delta s = l\cdot \Delta\theta=l\cdot\theta'\:\Delta t $.

La fuerza de fricción en el instante $t$ es $k\cdot \text{velocity} = k \cdot l\cdot \theta'(t)$. El trabajo total de fricción durante todo el intervalo es la integral del producto de la fuerza por el desplazamiento. $$ \text{Work} = \int \text{Friction }\: d s = \int_{t_0}^{t_f} \text{Friction }l \theta'\: d t= \int_{t_0}^{t_f} k l^2 (\theta')^2 \:d t $$ Aproxima la integral anterior para estimar la energía total disipada, usando algún método de cuadratura numérica. Explicita el método elegido, y explica el motivo de tu elección.

Pista: Antes de llamar al método de cuadratura, es conveniente obtener valores de $\theta'$ en una serie de puntos equiespaciados, y para ello puedes usar cualquier método de interpolación, pero también los atributos t_eval en python, o tspan en matlab.

  • En primer lugar resolvemos la ecuación con los valores dados, pero con el argumento tspan igualado a un linspace, para obtener los valores de $\theta$ y $\theta'$ en una serie de puntos equiespaciados.
In [16]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2); -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 30;
y0 = [0,1];

N = 2000;
tspan = linspace(t0, tf,N);

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

plot(ts, thetas, 'b-')
xlabel('t')
ylabel('theta')
title('Las oscilaciones se amortiguan...')
  • Ahora podemos usar cualquier método de integración que use un paso fijo: punto medio, trapecio, Simpson, etc...
In [17]:
function I = simpsons(f,a,b,n)
    % This function computes the integral "I" via Simpson's rule in the interval [a,b] with n+1 equally spaced points
    % 
    % Syntax: I = simpsons(f,a,b,n)
    % 
    % Where,
    %  f= can be either an anonymous function (e.g. f=@(x) sin(x)) or a vector
    %  containing equally spaced values of the function to be integrated
    %  a= Initial point of interval
    %  b= Last point of interval
    %  n= # of sub-intervals (panels), must be integer
    % 
    %  Written by Juan Camilo Medina  - The University of Notre Dame
    %  09/2010 (copyright Dr. Simpson)
    % 
    % 
    % Example 1:
    % 
    % Suppose you want to integrate a function f(x) in the interval [-1,1].
    % You also want 3 integration points (2 panels) evenly distributed through the
    % domain (you can select more point for better accuracy).
    % Thus:
    %
    % f=@(x) ((x-1).*x./2).*((x-1).*x./2);
    % I=simpsons(f,-1,1,2)
    % 
    % 
    % Example 2:
    % 
    % Suppose you want to integrate a function f(x) in the interval [-1,1].
    % You know some values of the function f(x) between the given interval,
    % those are fi= {1,0.518,0.230,0.078,0.014,0,0.006,0.014,0.014,0.006,0}
    % Thus:
    %
    % fi= [1 0.518 0.230 0.078 0.014 0 0.006 0.014 0.014 0.006 0];
    % I=simpsons(fi,-1,1,[])
    %
    % note that there is no need to provide the number of intervals (panels) "n",
    % since they are implicitly specified by the number of elements in the
    % vector fi
    if numel(f)>1 % If the input provided is a vector
        n=numel(f)-1;
        h=(b-a)/n;
        I= (h/3)*(f(1)+2*sum(f(3:2:end-2))+4*sum(f(2:2:end))+f(end));
    else % If the input provided is an anonymous function
        h=(b-a)/n; xi=a:h:b;
        I= (h/3)*(f(xi(1))+2*sum(f(xi(3:2:end-2)))+4*sum(f(xi(2:2:end)))+f(xi(end)));
    end
end
In [18]:
friction_force = k*l*vs;
work = friction_force.*vs;
%friction_work con metodo de integracion cutre (punto izquierdo)
'regla del punto izquierdo'
dt = (tf - t0)/N;
total_friction_work = dt*l*sum(work(1:N-1))
%friction_work con metodo del trapecio
'regla del trapecio'
total_friction_work = (dt/2)*l*(2*sum(work)-work(1)-work(N))
%friction_work con metodo de Simpson
'regla de Simpson'
total_friction_work = l*simpsons(work, t0, tf, N)

ans = regla del punto izquierdo
total_friction_work =  999.24
ans = regla del trapecio
total_friction_work =  996.24
ans = regla de Simpson
total_friction_work =  996.74
  • Curiosidad: El trabajo de fricción es exactamente igual a la diferencia de energía entre el instante inicial y el final, que se puede calcular sin hacer ninguna integral. $$k\cdot \text{velocity} = k \cdot l\cdot \theta'(t) = \left. \frac{1}{2}m l^2\theta'(t)^2 + m g l (1 - \cos(\theta(t))) \right|_{t=t_0}^{t=t_f} $$
In [19]:
theta0 = thetas(1);
v0 = vs(1);
thetaf = thetas(N);
vf = vs(N);

energia_0 = 0.5*m*(l^2)*(v0^2) + m*g*l*(1-cos(theta0));
energia_f = 0.5*m*(l^2)*(vf^2) + m*g*l*(1-cos(thetaf));
%El resultado es ligeramente diferente del obtenido integrando la
% energia disipada, a menos que impongamos una tolerancia muy baja 
% a solve_ivp
energia_0 - energia_f

ans =  997.56
  • Curiosidad: hacemos una gráfica con la evolución temporal de los tres tipos de energía: potencial, cinética y disipación. Observamos:
    • la disipación siempre aumenta.
    • la disipación aumenta más cuando la energía cinética es mayor.
    • la energía cinética tiene un máximo local cuando la potencial tiene un mínimo local, y viceversa.
    • la energía total es constante.
In [20]:
friction_work = dt*l*cumsum(friction_force.*vs);
kinetic_energy = 0.5*m*(l^2).*vs.^2;
potential_energy = m*g*l*(1-cos(thetas));
hold on;
plot(ts, kinetic_energy, 'b-')
plot(ts, potential_energy, 'g-')
plot(ts, friction_work, 'r-')
plot(ts, kinetic_energy + potential_energy + friction_work, 'k-')
hold off;
legend('kinetic', 'potential', 'losses due to friction', 'total energy, including friction')
xlabel('t')
ylabel('theta')
title('Energy breakdown by type...')
In [ ]: