Prescribing variables as an excitation in Runge-Kutta method

Question

I am using Runge-Kutta to solve a $3 \times 3$ 2^nd order linear ODE $$M x'' + C x' + K x = 0$$ and initial conditions are all zeros. Then I prescribe the 2nd variable to follow a given path.

As for prescribing displacement, I can choose two ways:

• [1] One is to use Lagrangian multiplier $$M x'' + C x' + K x + P x = P x_p,$$ where $x_p$ is prescribed path, it works perfectly.

• [2] Another way is to replace the variables with given values after Runge-Kutta updating, which will only work if the matrix $M$ is diagonal.

I really appreciate if anyone can offer any insights about this.

Even if I am prescribing only the displacement, it ends up that in Runge-Kutta, I have to prescribe the displacement rate as well, to match results with FD and other schemes.

I don't want to use the Lagrangian multiplier in Runge-Kutta or explicit methods, because it will significantly reduce my stability region and I have to use a very small time step.

Please feel free to jump with any suggestions.

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

My codes are attached below in case anyone needs to try. You can change $M$ matrix to compare the results.

Following is my main code. I put some parameters in a data structure so it's easier to access. And I have a prescribed displacement profile defined as 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);

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

Other scripts are given as well:

(Runge-Kutta Lagrangian multiplier)

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 - just replacing displacements after updating)

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

Below are the two model functions. The first one is the one used with the Lagrangian multiplier,

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

And the other model is used for replacing method,

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

Answer 1

In the friction-free case the Lagrangian formulation for the system with the given constraints would be $$ L(t,x,\dot x)=\frac12\dot x^TM\dot x-\frac12x^TKx+\lambda^T\Pi(x-x_p) $$ where the Lagrangian parameter $λ$ is also a function of $t$. $\Pi$ is the partial isometry that reduces the position to the components that get prescribed. The Euler Lagrange equations then give $$ 0=\frac{d}{dt}\frac{\partial L}{\partial \dot x}-\frac{\partial L}{\partial x} = M\ddot x+Kx-Π^Tλ. $$ The constraint $Πx=Πx_p$ leads to $Π\dot x=Π\dot x_p$ and $$ Π\ddot x_p=Π\ddot x=ΠM^{-1}(Π^Tλ-Kx)\impliesλ= (ΠM^{-1}Π^T)^{-1}Π(\ddot x_p+M^{-1}Kx) $$

I'm not sure what the "physical" way would be to incorporate the friction term.