# Prescribing variables as an excitation in Runge-Kutta method

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:

•  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.

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

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

• 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. – Lutz Lehmann May 5 '18 at 15:28
• 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. – Eric May 6 '18 at 13:58
• BTW, how do I @ someone when commenting? – Eric May 6 '18 at 14:00
• @Eric, you just add a @ before the username. – nicoguaro May 8 '18 at 16:12

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)$$