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

2023 7 27
0 个回答
更新时间：2023 7 28
8 次查看（30 天）

0 个投票

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

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Torsten
Torsten 2023-7-27
编辑：Torsten 2023-7-27
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
Yahor Staselka 2023-7-28
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
Torsten 2023-7-28
编辑：Torsten 2023-7-28
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
Yahor Staselka 2023-7-28
That's done the trick, thank you very much, Thorsten !!!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

产品

版本

R2021b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by