分数质量弹簧阻尼系统数学模型的数值解Matlab实现.rar

function[t,x,v]= FMSDS(m,b,k,alpha,beta,tau,dtau,x0,v0,f) % % DESCRIPTION: % % Function returns numerical solution of the % Fraction Mass-Spring-Damper System (FMSDS) % using Grunwald-Letnikov derivative definition. % % | % /| Damper -------- % /| -- | | % /|--- |-------| Mass | Force % /| -- b | m |<--------> f % /| Spring | | % /|-/\/\/\/\---| | % /| k -------- % /| | Displacement, x % /|<------------------------------> % /| | Velocity, v % 0 % % Fractional differential equation: % % d^alpha d^beta % m ---------- x(t) + b -------- x(t) + k x(t) = f(t) for 0 <= t <= tau % d t^alpha d^beta alpha > beta % % Initial conditions: % x(0) = x0 % x'(0) = v(0) = v0 % % Grunwald-Letnikov definition: % % Nf % sum b_j.y(t-j.dtau) % d^alpha j=0 % --------- y(t) = ------------------- % d t^alpha dtau^alpha % % Difference equation : % % p p % sum bc_j.x(p-j) sum cc_j.x(p-j) % j=0 j=0 % m --------------- + b --------------- + k.x(p) = f(p) % dtau^alpha dtau^beta % % or % p % m.dtau^(-alpha) ( bc_0.x(p) + sum bc_j.x(p-j) ) + % j=1 % % p % + b.dtau^(-beta) ( cc_0.x(p) + sum cc_j.x(p-j) ) + k.x(p) = f(p) % j=1 % % or % p p % x(p)=(f(p)-m.dtau^(-alpha).sum bc_j.x(p-j)-b.dtau^(-beta).sum cc_j.x(p-j))/ % j=1 j=1 % % (m.dtau^(-alpha).bc_0 + b.dtau^(-beta).cc_0 + k) % % INPUTS: % m - mass [ kg ] % b - damper coefficient [ kg/s ] % k - stiffness (or spring constant) [ kg/s^2 ] % alpha - fractional-order derivative [ - ] % beta - fractional-order derivative [ - ] % tau - simulation time [ s ] % dtau - time step [ s ] % x0 - initial displacement [ m ] % v0 - initial velocity [ m/s ] % f - force [ N ] % % OUTPUTS: % t - vector of time [ s ] % x - vector of distance [ m ] % v - vector of velocity [ m/s ] % % AUXILIARY VARIABLES: % n - length of time vector % p, i, j - index % bc, cc - binomial coefficients % s1, s2 - sum % % Copyright (C) 2020, Technical University of Kosice % Author : Terpak Jan % Revision: 17.02.2020 % n = tau/dtau+1; t(1) = 0.0; x(1) = x0; v(1) = v0; t(2) = dtau; x(2) = x(1) + dtau*v(1); v(2) = ( x(2) - x(1) )/dtau; p = 1; for i = 3:n p = p + 1; t(i) = (i-1)*dtau; bc = zeros(1,p+1); cc = zeros(1,p+1); bc(1) = 1; cc(1) = 1; for j = 2:p+1 bc(j) = (1 - (1+alpha)/(j-1))*bc(j-1); cc(j) = (1 - (1+beta)/(j-1))*cc(j-1); end s1 = 0; s2 = 0; for j = 2:p+1 s1 = s1 + bc(j)*x(p+2-j); s2 = s2 + cc(j)*x(p+2-j); end x(i) = ( f(i) - m*dtau^(-alpha)*s1 - b*dtau^(-beta)*s2 )/... ( m*dtau^(-alpha)*bc(1) + b*dtau^(-beta)*cc(1) + k); v(i) = ( x(i) - x(i-1) )/dtau; end end