在 Runge-Kutta 方法中规定变量作为激励

计算科学 matlab 颂 龙格库塔
2021-12-07 19:34:38

我正在使用 Runge-Kutta 解决一个3×3二阶线性ODE

Mx″+Cx′+Kx=0
初始条件全为零。然后我规定第二个变量遵循给定的路径。

至于规定位移，我可以选择两种方式：

  • [1] 一种是使用拉格朗日乘其中是规定路径，效果很好。

    Mx″+Cx′+Kx+Px=Pxp,
    xp

  • [2] 另一种方法是在 Runge-Kutta 更新后用给定值替换变量，这仅在矩阵是对角线时才有效。M

如果有人能对此提供任何见解，我真的很感激。

即使我只规定位移，最终在 Runge-Kutta 中，我也必须规定位移率，以使结果与 FD 和其他方案相匹配。

我不想在 Runge-Kutta 或显式方法中使用拉格朗日乘数，因为它会显着降低我的稳定性区域并且我必须使用非常小的时间步长。

请随时提出任何建议。

==================================================== ====================

我的代码附在下面，以防有人需要尝试。您可以更改矩阵来比较结果。M

以下是我的主要代码。我将一些参数放在数据结构中，以便更容易访问。我有一个规定的位移曲线定义为ExtDisp。

clc;
clear variables;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spring damper system

M = [3, 0.2 0.2;
    0.2, 2  0;
    0.2,    0,  5];

K = [10, 8 0.2;
    8, 10, 0.2;
    0.2, 0.2, 15];

C = zeros(3, 3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tspan
dt = 0.005;
t_span = [0, 8];
time = t_span(1): dt: t_span(2);

Ndt = length(time);
dim = length(M);

% which displacement to prescribe
index_prescribed = 2;
index_free = (1:dim)';
index_free(index_prescribed) = [];

% excitation path, a toneburst signal
freq = 1;
nCnt = 3.0;
ExtDisp_time_lim = nCnt / freq;
ExtDisp = zeros(1, Ndt);
for j = 1: Ndt
    if time(j) <= ExtDisp_time_lim
        ExtDisp(j) = (1-cos(2*pi*time(j)*freq/nCnt)) * sin(2*pi*time(j)*freq);
    end
end

dExtDisp = diff(ExtDisp) / dt;
ddExtDisp = diff(dExtDisp) / dt;

force_param_cells.time = time;
force_param_cells.ExtDisp_lim = ExtDisp_time_lim;
force_param_cells.ExtDisp = ExtDisp;
force_param_cells.dExtDisp = dExtDisp;
force_param_cells.ddExtDisp = ddExtDisp;
force_param_cells.index_prescribed = index_prescribed;

% MAIN COMPUTATIONS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use Lagrangian multiplier

% initial BCs
x0 = [0, 0, 0]';
v0 = [0, 0, 0]';
[tt2, yy2] = rk_lagrangian_multiplier(dt, t_span, M, C, K, force_param_cells, x0, v0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% just by replacing

% initial BCs
x0 = [0, 0, 0]';
v0 = [0, 0, 0]';
[tt4, yy4] = rk_ordinary(dt, t_span, M, C, K, force_param_cells, x0, v0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot all variables
figure;hold on;
for k = 1 : 6
    subplot(2, 3, k);hold on;
    plot(tt2, yy2(k,:), 'k-', 'LineWidth', 2);
    plot(tt4, yy4(k,:), 'r--', 'LineWidth', 2);

    if k > 3
        ylabel(['$\dot{x}$', int2str(k-3)], 'FontSize', 15 ,'interpreter','latex');
    else
        ylabel(['${x}$', int2str(k)], 'FontSize', 15 ,'interpreter','latex');
    end

end
legend('RK-Multiplier', 'RK-Replace');
xlabel('Time (s)','FontSize',15);

==================================================== =================

还给出了其他脚本：

（龙格-库塔拉格朗日乘数）

function [tt, yy] = rk_lagrangian_multiplier(dt, t_span, M, C, K, force_param_cells, x0, v0)

time = t_span(1): dt: t_span(2);

% system dimensions
dim = length(M);
Ndt = length(time);

index_prescribed = force_param_cells.index_prescribed;

% penalty on prescribed displacement
lamda = max(max(K)) * 1E4;
PenaltyMatrix = zeros(3);
PenaltyMatrix(index_prescribed, index_prescribed) = lamda;

% save model parameters as cells
model_cells{1} = M;
model_cells{2} = C;
model_cells{3} = K;
model_cells{4} = PenaltyMatrix;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Runge-Kutta %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

step = 1;
tt(:,step) = time(step);
yy(:,step) = [x0; v0];

for step = 1 : Ndt-1          % calculation loop
    time_t = time(step);
    yy_t = yy(:, step);

    k_1 = model_rk_spring_mass_full(time_t,              yy_t,                   model_cells,    force_param_cells);
    k_2 = model_rk_spring_mass_full(time_t + 0.5 * dt,  yy_t + 0.5*dt * k_1,    model_cells,    force_param_cells);
    k_3 = model_rk_spring_mass_full(time_t + 0.5 * dt,  yy_t + 0.5*dt * k_2,    model_cells,    force_param_cells);
    k_4 = model_rk_spring_mass_full(time_t + dt,         yy_t + dt * k_3,        model_cells,   force_param_cells);

    tt(step + 1)    =   time_t + dt;
    yy(:, step + 1) =   yy(:, step) + dt * (k_1 + 2*k_2 + 2*k_3 + k_4) / 6.0;
end

（Runge-Kutta - 更新后只替换位移）

function [tt, yy] = rk_ordinary(dt, t_span, M, C, K, force_param_cells, x0, v0)

time = t_span(1): dt: t_span(2);

% system dimensions
dim = length(M);
Ndt = length(time);

% save model parameters as cells
model_cells{1} = M;
model_cells{2} = C;
model_cells{3} = K;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Runge-Kutta %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

step = 1;
tt(:,step) = time(step);
yy(:,step) = [x0; v0];
yy(force_param_cells.index_prescribed, step) = force_param_cells.ExtDisp(:, step);

for step = 1 : Ndt-1          % calculation loop

    time_t = time(step);
    yy_t = yy(:, step);

    k_1 = model_rk_spring_mass_no_penalty(time_t,               yy_t,                   model_cells);
    k_2 = model_rk_spring_mass_no_penalty(time_t + 0.5 * dt,    yy_t + 0.5*dt * k_1,    model_cells);
    k_3 = model_rk_spring_mass_no_penalty(time_t + 0.5 * dt,    yy_t + 0.5*dt * k_2,    model_cells);
    k_4 = model_rk_spring_mass_no_penalty(time_t + dt,          yy_t + dt * k_3,        model_cells);

    tt(step + 1)    =   time_t + dt;
    yy(:, step + 1) =   yy(:, step) + dt * (k_1 + 2*k_2 + 2*k_3 + k_4) / 6.0;

    % prescribing displacement and rate
    if time_t <= force_param_cells.ExtDisp_lim
        yy(force_param_cells.index_prescribed, step+1) = force_param_cells.ExtDisp(step+1);
        yy(force_param_cells.index_prescribed + dim, step+1) = force_param_cells.dExtDisp(step+1);
    end

end

下面是两个模型函数。第一个是与拉格朗日乘数一起使用的，

function dy = model_rk_spring_mass_full(t, y, model_cells, force_param_cells)

M = model_cells{1};
C = model_cells{2};
K = model_cells{3};
Penalty = model_cells{4};

time = force_param_cells.time;
ExtDisp = force_param_cells.ExtDisp;
ExtDisp_lim = force_param_cells.ExtDisp_lim;
index_prescribed = force_param_cells.index_prescribed;

dim = length(M);
dy = zeros(dim * 2, 1);

% dealing with extra prescribed displacement

disp_prescribed = interp1(time', ExtDisp', t)';
disp_vec = zeros(dim, 1);
disp_vec(index_prescribed) = disp_prescribed;

if t <= ExtDisp_lim
    P_m = Penalty;
    force_penalty = P_m * disp_vec;
else
    P_m = zeros(dim, dim);
    force_penalty = zeros(dim, 1);
end

% updating procedures
dy(1:dim)       = y(dim+1:end);
dy(dim+1:end)   = M \ (force_penalty - C * y(dim+1: end) - (K + P_m) * y(1:dim));

end

而另一个模型用于替换方法，

function dy = model_rk_spring_mass_no_penalty(t, y, model_cells)

M = model_cells{1};
C = model_cells{2};
K = model_cells{3};

dim = length(M);
dy = zeros(dim * 2, 1);

% updating procedures
dy(1:dim)       = y(dim+1:end);
dy(dim+1:end)   = M \ (- C * y(dim+1: end) - K * y(1:dim));

end
1个回答

在无摩擦情况下，具有给定约束的系统的拉格朗日公式将是 其中拉格朗日参数也是的函数。是部分等距，它将位置减少到指定的组件。欧拉拉格朗日方程则给出 约束导致和

L(t,x,x˙)=12x˙TMx˙−12xTKx+λTΠ(x−xp)
λtΠ
0=ddt∂L∂x˙−∂L∂x=Mx¨+Kx−ΠTλ.
Πx=ΠxpΠx˙=Πx˙p
Πx¨p=Πx¨=ΠM−1(ΠTλ−Kx)⟹λ=(ΠM−1ΠT)−1Π(x¨p+M−1Kx)

我不确定包含摩擦术语的“物理”方式是什么。

其它你可能感兴趣的问题
上一篇如何为DFT制作部分占用的超级单元？ 下一篇收敛分析期间斜率减小