function dotq = beam_supported_dynamic_equation_param (t,q, n,M,B,K,constraints,alpha,force)
% ビームの運動方程式
% a:制約を表す行列
    un = q(1:n);
    vn = q(n+1:2*n);
    
    [am,an] = size(constraints);
    
    mat = [ M,  -constraints; ...
           -constraints',  zeros(an,an) ];
    vec = [ -K*un-B*vn + force(t,n,un,vn); ...
            constraints'*(2*alpha*vn + (alpha^2)*un) ];
    sol = mat \ vec;
    dotvn = sol(1:n);
    
    dotq = [ vn; dotvn ];
end
