0
$\begingroup$

I am using Runge-Kutta to solve a $3 \times 3$ 2nd 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
$\endgroup$
  • $\begingroup$ What exactly are you hoping to achieve? Make $x$ follow closely the path $x_p$ while still solving the ODE? Or have the free dynamic given by the first system bound to some surface? Note that the second case makes your ODE into a index 2 or index 3 DAE system so that indeed projection will violate the hidden constraints. Most simple DAE solvers are based on backwards-differentiation formula (BDF) multistep methods, higher index systems require index analysis and index reduction to semi-explicit index 1 or 2 systems. $\endgroup$ – LutzL May 5 '18 at 15:28
  • $\begingroup$ Well, thanks for pointing this out. I am solving a system with one variable following a path, which is quite common in some physics. Yes, reducing the system dimension will work, and it will match the results with using Lagrangian multiplier. But replacing or projection (as you mentioned) will work in some scenarios as well, when the matrix M is diagonal. $\endgroup$ – Eric May 6 '18 at 13:58
  • $\begingroup$ BTW, how do I @ someone when commenting? $\endgroup$ – Eric May 6 '18 at 14:00
  • 1
    $\begingroup$ @Eric, you just add a @ before the username. $\endgroup$ – nicoguaro May 8 '18 at 16:12
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Hi Lutzl, thanks for commenting on my questions. What you presented basically talks about how to address constraints for a physical system. But if you consider it mathematically, it doesn't matter whether there is friction term or not, you just give Πx=Πxp, on the left and right ends of your equation. But note that Π needs to be several orders larger than system matrices. $\endgroup$ – Eric May 7 '18 at 20:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.