布料仿真先导3-带阻尼的双球单摆下的拉格朗日方程列些和matlab仿真

这个是在先导2基础上，将单小球系统变为双小球系统，每个小球都会受到各自的阻尼力而运动。

 

matlab代码如下

%基于ode45的双球摆动带阻尼的运动方程
%此处是基于拉个朗日函数分析法给出的微分方程求解

close all;clearvars;clc;
dbstop if error;

global m1 m2 l1 l2 g k1 k2;
m1 = 1;
m2 = 1;
l1 = 1;
l2 = 1;
g = 9.8;
k1 = 1.2;
k2 = 1.5;

[t , y] = ode45(@vdp1,[0 20],[50*pi/180 0 40*pi/180 0]);

x1 = y(: , 1) .* 180/pi;
x2 = y(: , 3) .* 180/pi;
h1 = figure;
for i = 1 : 1 : length(t)
    set(0,'CurrentFigure',h1);
    scatter(0 , 0 , '*');hold on;
    x = [l1*sind(x1(i)) , l1*sind(x1(i)) + l2*sind(x2(i))];
    y = [-1*l1*cosd(x1(i)) , -l1*cosd(x1(i)) - l2*cosd(x2(i))];
    scatter(x , y);hold on;
    plot([0 x(1)] , [0 y(1)] , 'r');hold on;
    plot([x(1) x(2)] , [y(1) y(2)] , 'r');hold on;
    xlim([-3 3]);
    ylim([-3 3]);
    pause(0.05);
    hold off;
end

function dydt = vdp1(t,y)   %y = [phi1 phi1']输入为弧度
    
    global m1 m2 l1 l2 g k1 k2;
    dydt = zeros(4 , 1);
    dydt(1) = y(2);

    A = (m1 + m2)*l1;
    B = m2*l2*cos(y(1) - y(3));
    C = m2*l2*sin(y(1) - y(3));
    D = k2*l1*l1*cos(y(1) - y(3));
    E = (k1 + k2)*l1*l1;
    F = (m1+m2)*g*sin(y(1));
    H = m2*l2;
    I = m2*l1*cos(y(1) - y(3));
    J = m2*l1*sin(y(1) - y(3));
    K = k2*l2*l2;
    L = k2*l1*l2*cos(y(1) - y(3));
    M = m2*g*sin(y(3));
    dydt(2) = ((B/H)*(J*y(2)^2+K*y(4)+L*y(2)+M) - (C*y(4)*y(4)+D*y(4)+E*y(2)+F))/(A-B*I/H);
    dydt(3) = y(4);
    dydt(4) = ((-1*M-L*y(2)-K*y(4)-J*y(2)*y(2)) + (I/A)*(F+E*y(2)+D*y(4)+C*y(4)*y(4)))/(H-I*B/A);
 
end