微分方程数值算法求解结构动力响应

八达鸟007

已于 2024-07-13 19:51:36 修改

阅读量401

点赞数 3

分类专栏：龙格-库塔方法欧拉法文章标签：算法 matlab

于 2024-07-13 17:12:53 首次发布

本文链接：https://blog.csdn.net/qq_33028799/article/details/140400382

版权

龙格-库塔方法同时被 2 个专栏收录

1 篇文章 0 订阅

订阅专栏

欧拉法

1 篇文章 0 订阅

订阅专栏

对于如下的二自由度结构，采用不同的微分方程解法来求解结构动力响应。

在这里插入图片描述
通过牛顿第二定律，可以得到如下的动力学方程：
$\left[ \begin{matrix} m_1& 0\\ 0& m_1\\ \end{matrix} \right] \left\{ \begin{array}{c} \ddot{x}_1\\ \ddot{x}_2\\ \end{array} \right\} +\left[ \begin{matrix} c_1+c_2& -c_2\\ -c_2& c_2\\ \end{matrix} \right] \left\{ \begin{array}{c} \dot{x}_1\\ \dot{x}_2\\ \end{array} \right\} +\left[ \begin{matrix} k_1+k_2& -k_2\\ -k_2& k_2\\ \end{matrix} \right] \left\{ \begin{array}{c} x_1\\ x_2\\ \end{array} \right\} =\left\{ \begin{array}{c} f_1\\ f_2\\ \end{array} \right\}$
简化为：
$\mathbf{M\ddot{X}}+\mathbf{C\dot{X}}+\mathbf{KX}=\mathbf{F}$
其中：
$\mathbf{M}=\left[ \begin{matrix} m_1& 0\\ 0& m_1\\ \end{matrix} \right] ,\mathbf{C}=\left[ \begin{matrix} c_1+c_2& -c_2\\ -c_2& c_2\\ \end{matrix} \right] ,\mathbf{K}=\left[ \begin{matrix} k_1+k_2& -k_2\\ -k_2& k_2\\ \end{matrix} \right] , \mathbf{X}=\left\{ \begin{array}{c} x_1\\ \,\,x_2\\ \end{array} \right\}$
写为状态空间方程，可得
$\left\{ \begin{array}{c} \mathbf{\dot{X}}\\ \mathbf{\ddot{X}}\\ \end{array} \right\} =\left[ \begin{matrix} \mathbf{0}& \mathbf{I}\\ -\mathbf{M}^{-1}\mathbf{C}& -\mathbf{M}^{-1}\mathbf{K}\\ \end{matrix} \right] \left\{ \begin{array}{c} \mathbf{X}\\ \mathbf{\ddot{X}}\\ \end{array} \right\} +\left[ \begin{array}{c} \bf{0}\\ \mathbf{M}^{-1}\\ \end{array} \right] \mathbf{F}$

简化为：
$\mathbf{\dot{Z}}=\mathbf{AZ}+\mathbf{BU}$
其中：

$\mathbf{Z}=\left\{ \begin{array}{c} \mathbf{X}\\ \mathbf{\dot{X}}\\ \end{array} \right\} \text{，}\mathbf{A}=\left[ \begin{matrix} \mathbf{0}& \mathbf{I}\\ -\mathbf{M}^{-1}\mathbf{C}& -\mathbf{M}^{-1}\mathbf{K}\\ \end{matrix} \right] \text{，}\mathbf{B}=\left[ \begin{array}{c} \mathbf{0}\\ \mathbf{M}^{-1}\\ \end{array} \right] \text{，}\mathbf{U}=\mathbf{F}$

上式为一微分方程，可采用微分多种数值解法来进行求解，包括欧拉法、改进欧拉法、龙格-库塔法、亚当姆斯法、亚当姆斯预报-校正法等，以下是相应的MATLAB代码.


function Y = state_space_slover(A,B,C,D,U,t,Z0,method)


switch method
    case 1
        Y = Euler_method(A,B,C,D,Z0,U,t);
    case 2
        Y = Improved_Euler_method(A,B,C,D,Z0,U,t);
    case 3
        Y = Runge_kutta_method_4th(A,B,C,D,Z0,U,t);
    case 4
        Y = Adams_method_2nd(A,B,C,D,Z0,U,t);
    case 5
        Y = Adams_method_4th(A,B,C,D,Z0,U,t);
    case 6
        Y = Improver_Adams_method_4th(A,B,C,D,Z0,U,t);
end

end

function Y = Euler_method(A,B,C,D,Z0,U,t)
dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    z0 = Z(:,it);
    Ui = U(:,it);
    zi = A*z0 + B*Ui;
    
    Z(:,it+1) = z0 + dt*zi;
    Y(:,it) = C*zi + D*Ui;
end

end

function Y = Improved_Euler_method(A,B,C,D,Z0,U,t)

dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    z0 = Z(:,it);
    Ui = U(:,it);

    zp = z0 + dt*(A*z0 + B*Ui);
    zc = z0 + dt*(A*zp + B*Ui);
    zi = (zp + zc)/2;

    Z(:,it+1) = zi;
    Y(:,it) = C*zi + D*Ui;
end

end

function Y = Runge_kutta_method_4th(A,B,C,D,Z0,U,t)

dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    z0 = Z(:,it);
    Ui = U(:,it);

    K1 = A*z0 + B*Ui;
    K2 = A*(z0 + K1*dt/2) + B*Ui;
    K3 = A*(z0 + K2*dt/2) + B*Ui;
    K4 = A*(z0 + K3*dt) + B*Ui;
    zi = z0 + dt/6*(K1 + 2*K2 + 2*K3 + K4);

    Z(:,it+1) = zi;
    Y(:,it) = C*zi + D*Ui;
end

end

function Y = Adams_method_2nd(A,B,C,D,Z0,U,t)

dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    Ui = U(:,it);
    z0 = Z(:,it);
    if it == 1
        zi = A*Z(:,it) + B*Ui;
    else
        z1 = Z(:,it);
        z2 = Z(:,it-1);

        za = A*z1 + B*Ui;
        zb = A*z2 + B*Ui;
        zi = z0  + 0.5*dt*(3*za - zb);
    end

    Z(:,it+1) = zi;
    Y(:,it) = C*zi + D*Ui;

end

end

function Y = Adams_method_4th(A,B,C,D,Z0,U,t)

dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    Ui = U(:,it);
    z0 = Z(:,it);
    if it <= 4
        K1 = A*z0 + B*Ui;
        K2 = A*(z0 + K1*dt/2) + B*Ui;
        K3 = A*(z0 + K2*dt/2) + B*Ui;
        K4 = A*(z0 + K3*dt) + B*Ui;
        zi = z0 + dt/6*(K1 + 2*K2 + 2*K3 + K4);
    else
        z1 = Z(:,it);
        z2 = Z(:,it-1);
        z3 = Z(:,it-2);
        z4 = Z(:,it-3);

        za = A*z1 + B*Ui;
        zb = A*z2 + B*Ui;
        zc = A*z3 + B*Ui;
        zd = A*z4 + B*Ui;

        zi = z0  + 1/24*dt*(9*za + 19*zb - 5*zc + zd);
    end
    Z(:,it+1) = zi;
    Y(:,it) = C*zi + D*Ui;
end
end

function Y = Improver_Adams_method_4th(A,B,C,D,Z0,U,t)

dt = t(2) - t(1);
Nt = length(t);
Z = zeros(size(A,1), Nt);
Z(:,1) = Z0;
Y = zeros(size(C,1), Nt);

for it = 1:Nt
    Ui = U(:,it);
    z0 = Z(:,it);
    if it <= 3
        K1 = A*z0 + B*Ui;
        K2 = A*(z0 + K1*dt/2) + B*Ui;
        K3 = A*(z0 + K2*dt/2) + B*Ui;
        K4 = A*(z0 + K3*dt) + B*Ui;
        zi = z0 + dt/6*(K1 + 2*K2 + 2*K3 + K4);
    else
        z1 = Z(:,it);
        z2 = Z(:,it-1);
        z3 = Z(:,it-2);
        z4 = Z(:,it-3);

        za = A*z1 + B*Ui;
        zb = A*z2 + B*Ui;
        zc = A*z3 + B*Ui;
        zd = A*z4 + B*Ui;

        ze = z0 + dt/24*(55*za - 59*zb + 37*zc - 9*zd);
        zp = A*ze + B*Ui;

        zi = z0 + dt/24*(9*zp + 19*za - 5*zb + zc);
    end

    Z(:,it+1) = zi;
    Y(:,it) = C*zi + D*Ui;

end

end

假设 $m_1$ = 10 kg, $m_2$ = 5kg, $c_1$ = 3 $\text{N}\cdot \text{s}/\text{m}$ , $c_2= 4 \text{N}\cdot \text{s}/\text{m}$ , , $k_1 = 10 \text{N}/\text{m}$ , $k_2$ = 6 $\text{N}/\text{m}$ , $f_1 = 5\sin(t)$ , $f_2 = 10\sin(t)$ 。若只考虑输出结构的位移和速度响应，则状态空间方程的数值为：

$\mathbf{A}=\left[ \begin{matrix} 0& 0& 1& 0\\ 0& 0& 0& 1\\ -1.6& 0.6& -0.7& 0.4\\ 1.2& -1.2& 0.8& -0.8\\ \end{matrix} \right] ,\mathbf{B}=\left[ \begin{matrix} 0& 0\\ 0& 0\\ 0.1& 0\\ 0& 0.2\\ \end{matrix} \right] ,\mathbf{C}=\left[ \begin{matrix} 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ \end{matrix} \right] ,\mathbf{D}=\left[ \begin{matrix} 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ \end{matrix} \right]$

由此，采用以上方法，可得到结构的响应值，代码如下：


clear
clc

m1 = 10;
m2 = 5;
c1 = 3;
c2 = 4;
k1 = 10;
k2 = 6;

M = [m1, 0; 0, m2];
C = [c1+c2, -c2; -c2, c2];
K = [k1+k2, -k2; -k2, k2];

n = 2;
A = [zeros(n), eye(2); -M\K, -M\C];
B = [zeros(n); inv(M)];
C = eye(2*n);
D = zeros(2*n,n);

t = 0:0.01:10;
f1 = 5*sin(t);
f2 = 10*sin(t);
U = [f1; f2];
Z0 = [1,0.5,0,0]';

Y = cell(1,6);
for i = 1:6
    Y{i} = state_space_slover(A,B,C,D,U,t,Z0,i);
end

于此同时，可采用MATLAB中 lsim 命令求得此状态方程的精确值

Y0 = lsim(A,B,C,D,U,t,Z0);

不同计算方式的结果对比如下：

colors = {'#a9eee6', '#ffbd39', '#f38181', '#fccde2', '#5e63b6', '#685454'};
style = {'-','-.', '--', '.', '-.', '-.'};
ylabels = {'$x_1$', '$x_2$', '$\dot{x}_1$', '$\dot{x}_2$'};
figure(Color='w', Position=[200,200,800,600])
for i = 1:4
    subplot(2,2,i)

    plot(t,Y0(:,i),'b',linewidth = 1.5)
    hold on
    for j = 1:6
        plot(t, Y{j}(i,:), Color=colors{j}, linewidth = 1.5)
    end

    if i == 1
        legend('lsim', 'Euler', 'Improved Euler', '4th Runge kutta', '2nd Adams', ...
            '4th Adams','4th Improved Adams', ...
            'Location', 'NorthOutside', 'Orientation', 'horizontal','NumColumns',4) 
    end

    xlabel('Time (s)')
    ylabel(ylabels{i}, 'Interpreter','latex',FontSize=18)
end