Solving ODEs in a for-loop with an additional variable.

Hi!

I'm trying to solve a system of 2 ODEs in a for-loop. Each iteration of the for-loop changes a variable in the ODE. My issue is that the ode45 is taking only the first value of this variable and discards the rest.

What is the correct way of iterating through parameters "xr" and "vr" in the last function?

Sorry for a newbie question, any help will be greatly appreciated!

xr = amplitude*sin(2*pi*time*heave_input_frequency); % displacement

vr = diff(xr)/SIM_SETUP.time_step; % speed

ar = diff(vr)/SIM_SETUP.time_step; % acceleration

output(:,2)=xr; % Record displacement in the 2nd column of the output table

output(1:end-1,3)=vr;output(end,3)=output(end-1,3); % Record speed in the 3rd column of the output table

output(1:end-2,4)=ar;output(end-2:end,4)=output(end-2,4); % Record speed in the 4th column of the output table

for i=1:length(time)

% -------------------- Retrieve input

xr=output(i,2);

vr=output(i,3);

[t,q] = ode45(@(t,q) f(t,q,CAR_SETUP,xr,vr),time,y0,opts);

% -------------------- Record data

output(i,5)=q(i,1); % xs

output(i,6)=q(i,2); % vs

output(i,7)=q(i,3); % xt

output(i,8)=q(i,4); % vt

% -------------------- Next iteration

i=i+1;

end

function M = mass(t,q,CAR_SETUP,xr,vr)

% Extract parameters

ms = CAR_SETUP.ms;

mt = CAR_SETUP.mt;

% Mass matrix function

M = zeros(4,4);

M(1,1) = 1;

M(2,2) = ms;

M(3,3) = 1;

M(4,4) = mt;

end

function dydt = f(t,q,CAR_SETUP,xr,vr)

% Extract parameters

ks = CAR_SETUP.ks;

cs = CAR_SETUP.cs;

kt = CAR_SETUP.kt;

ct = CAR_SETUP.ct;

% Equation to solve

dydt = [q(2)

cs*(q(4)-q(2))+ks*(q(3)-q(1))

q(4)

-cs*(q(4)-q(2))+ct*(vr-q(4))-ks*(q(3)-q(1))+kt*(xr-q(3))];

end

[Torsten]

We cannot run your code because several variables are undefined. Thus we are not able to answer your question. Especially setting i = i +1 and writing a single value of the solution vector into the output array seems strange:

output(i,5)=q(i,1);                             % xs
output(i,6)=q(i,2);                             % vs
output(i,7)=q(i,3);                             % xt
output(i,8)=q(i,4);                             % vt

% -------------------- Next iteration
i=i+1;

[Yahor Staselka]

Apologies, Thorsten, here is the full code. In the for loop i'm only writing the results of the ODE solution for that particular iteration of the loop, hence the single values.

I'm guessing my whole approach is not correct in this case. I have seen systems of ODEs solved this way with only the initial conditions, but no additional inputs. In my case these inputs are different for every iteration of the loop and i'm not quite sure how to include them...

I've tried adding them as additional parameters, but that only seemed to work for me when these values are constant, not changing every step of the iteration.

Many thanks for your time!

clear all;close all;clc;tic;

%% -------------------- Simulation setup
SETUP.ms = 347.2;
SETUP.ks = 40000;
SETUP.cs = 4500;
SETUP.mt = 52.8;
SETUP.kt = 312000;
SETUP.ct = 5;

calc_frequency = 100;           % Hz
time_step = 1/calc_frequency;   % s
sim_time = 2;                   % s

opts = odeset('Mass',@(t,q) mass(t,q,SETUP));
time=(0:time_step:sim_time);

%% -------------------- Prepare output table
output=zeros(length(time),8);   % Empty output table
output(:,1)=time;               % Record time in the 1st column of the output table

%% -------------------- Generate road input
amplitude = 0.1;                                                % m, amplitude
heave_input_frequency = 2;                                      % Hz
xr = amplitude*sin(2*pi*time*heave_input_frequency);            % input displacement
vr = diff(xr)/time_step;                                        % input speed
ar = diff(vr)/time_step;                                        % input acceleration
output(:,2)=xr;                                                 % Record input displacement in the 2nd column of the output table
output(1:end-1,3)=vr;output(end,3)=output(end-1,3);             % Record input speed in the 3rd column of the output table
output(1:end-2,4)=ar;output(end-2:end,4)=output(end-2,4);       % Record input speed in the 4th column of the output table

%% -------------------- Simulation
i=1;
y0 = [0; 0; 0; 0];
progress_bar = waitbar(0, 'Starting');

for i=1:length(time)

    % -------------------- Retrieve input
    xr=output(i,2);
    vr=output(i,3);

    [t,q] = ode45(@(t,q) f(t,q,SETUP,xr,vr),time,y0,opts);

    % -------------------- Record data
    output(i,5)=q(i,1);                             % xs
    output(i,6)=q(i,2);                             % vs
    output(i,7)=q(i,3);                             % xt
    output(i,8)=q(i,4);                             % vt
   
    % -------------------- Progress bar update
    waitbar(i*time_step/sim_time, progress_bar, ...
        sprintf('Progress: %d %%', floor(i*time_step/sim_time*100)));
    
    % -------------------- Next iteration
    i=i+1;
end

%% -------------------- Plot the results
t =output(:,1);
xr=output(:,2);
vr=output(:,3);
ar=output(:,4);
xs=output(:,5);
vs=output(:,6);
xt=output(:,7);
vt=output(:,8);

subplot(3,1,1)
hold on
plot(t,xr);
plot(t,xs);
plot(t,xt);
legend('Disp road','Disp chassis','Disp tyre')
hold off
grid on
xlabel('Time (s)')
ylabel('Displacement (m)')

toc

%% -------------------- Functions
function M = mass(t,q,SETUP,xr,vr)
    % Extract parameters
    ms = SETUP.ms;
    mt = SETUP.mt;
    % Mass matrix function
    M = zeros(4,4);
    M(1,1) = 1;
    M(2,2) = ms;
    M(3,3) = 1;
    M(4,4) = mt;
end

function dydt = f(t,q,SETUP,xr,vr)
    % Extract parameters
ks = SETUP.ks;
cs = SETUP.cs;
kt = SETUP.kt;
ct = SETUP.ct;
    
% Equation to solve
dydt = [q(2)
        cs*(q(4)-q(2))+ks*(q(3)-q(1))
        q(4)
        -cs*(q(4)-q(2))+ct*(vr-q(4))-ks*(q(3)-q(1))+kt*(xr-q(3))];
end

[Torsten]

Generate three numerical arrays: One for a time vector (T) and the other two for xr and vr (XR and VR), evaluated at these times specified in the vector T.

Then define interpolation functions

xr = @(t) interp1(T,XR,t)
vr = @(t) interp1(T,VR,t)

pass these functions to f and evaluate them at the time t that is input from the solver:

function dydt = f(t,q,SETUP,xr,vr)
% Extract parameters
ks = SETUP.ks;
cs = SETUP.cs;
kt = SETUP.kt;
ct = SETUP.ct;

% Equation to solve
dydt = [q(2)
        cs*(q(4)-q(2))+ks*(q(3)-q(1))
        q(4)
        -cs*(q(4)-q(2))+ct*(vr(t)-q(4))-ks*(q(3)-q(1))+kt*(xr(t)-q(3))];
end

If you have further problems, look at the example

ODE with Time-dependent Terms

under

https://uk.mathworks.com/help/matlab/ref/ode45.html

[Yahor Staselka]

That's done the trick, thank you very much, Thorsten !!!