How to solve this ode equation where one of the inputs has some delay and equation has written in matrix form in MATLAB ode45?

I solved this ode for constant input u. Function script of this ode is written here-

function dx = mySystem(t,x,u,A,B)
dx = A*x +B*u
end

A is (6*6) matrix. B is (4*6) matrix. C is (6*1) matrix. D is (4*1) matrix.

Main Script

A=[ -0.1078 0 0 0 -10 0; 0 -0.0898 0 0 0 -8.3333; 0 0 -0.2352 0 10 0; 0 0 0 -0.196 0 8.3333; 1 0 -1 0 0 0; 0 1 0 -1 0 0];
B=[0.3313 0 -8.9335 0; 0 0.2761 0 -7.4446; 0 0 3.7249 0;0 0 0 3.104; 0 0 0 0; 0 0 0 0];
C=[ 0 0 1.4206 1.4206 0 0];  D= [0 0 -15 -15];

u=zeros(2*2,1);
for i=1:2
    u(i)=0.5;
end
fun = @(t, x) mySystem(t, x, u, A, B);
tspan = [0, 20]; % Time span
x0 = zeros(3*3,1); %Initial Condition

[t, YSol] = ode45(fun, tspan, x0);  % ode

Pe=C*YSol(:,1:3*2)'+D*u; % Output Equation
% Plot
figure(5)
plot(t,Pe);
xlabel('time'); ylabel('{P');
title(' Change in P')

Above ode has been solved easily. But actually my system has delayed input. u(1)=0.5 is constant but u(2)=0.5*(tau-5), i.e. u(2) is delayed by 10 second. In this situation how can I solve this? What will be the modified code? Other thing is if the system have more than 100 order then how it can be solved? A is (3*n by3*n); B is(3*n by 2*n); C is (1 by 3*n);D is (1 by 2*n); where n=100, and out of 2*n input u(2):u(n) have delay of different tau.

%% Delay in Input 
u(1)=0.5;

if t <= 10
    u(2)=0;
else 
     u(2)=0.5;
end

[Walter Roberson]

See the discussion at https://www.mathworks.com/matlabcentral/answers/1671864-trying-to-use-ode45-to-solve-a-second-order-differential-equation-with-time-dependent-parameters#answer_917859

Answer 1

[Torsten]

Order the complete list of time instances when the control variables change their values.

If these times are 0 < t1 < t2 < t3 < ... < tn, then call ode45 in a loop and integrate first with tspan1 = [0,t1], then with tspan2 = [t1,t2] by using the end values at t1 as start values at t2 and so on.

[Sol Elec]

Thanks a lot @Torsten. Here is the my different time ode equation. How to form a loop for this ode solution?

% Simulation run time 0 to 40 seconds
Delay = [0; 3; 6; 9; 12]; 
tend=25;
vpu =[0.0417;0.0394;0.0393;0.0392;0.0392];

u(1)=vpu(1);
tspan1 = [Delay(1), Delay(2)]; % Time Span
x0 = zeros(3*n,1); %Initial Condition
[t1, YSol1] = ode45(@(t, x) myFun(t, x, u, A, B), tspan1, x0);

u(2)=vpu(2);
tspan2 = [Delay(2), Delay(3)];
[t2, YSol2] = ode45(@(t, x) myFun(t, x, u, A, B), tspan2, YSol1(end, :));

u(3)=vpu(3);
tspan3 = [Delay(3), Delay(4)];
[t3, YSol3] = ode45(@(t, x) myFun(t, x, u, A, B), tspan3, YSol2(end, :));

u(4)=vpu(4);
tspan4 = [Delay(4), Delay(5)];
[t4, YSol4] = ode45(@(t, x) myFun(t, x, u, A, B), tspan4, YSol3(end, :));

u(5)=vpu(5);
tspan5 = [Delay(5), tend];
[t5, YSol5] = ode45(@(t, x) myFun(t, x, u, A, B), tspan5, YSol4(end, :));

t=[t1;t2;t3;t4;t5];
YSol=[YSol1;YSol2;YSol3;YSol4;YSol5];

[Torsten]

The code will work out fine as is.

Do you need it in a loop because your real problem has many more such points of interruption (i.e. the delay vector will be much longer) ?

[Sol Elec]

Yes I need it in loop, but how to convert it into loop? The main problen is related to t in loop [t1, YSol1].

[Torsten]

tend = 25;
Delay = [0 ;3; 6; 9; 12; tend]; 
vpu =[0.0417;0.0394;0.0393;0.0392;0.0392];
x0 = zeros(3*n,1);   % initialize initial conditions
T = [];
X = [];
for i = 1 : numel(vpu)
    u = vpu(i);
    tspan = [Delay(i) , Delay(i+1)]; % Time Span
    [t , x] = ode45(@(t,x) myfun(t, x, u, A, B), tspan, x0);
    x0 = x(end , :);
    T = [T ; t(1:end-1)];
    X = [X ; x(1:end-1,:)];
end

[Sol Elec]

Thanks a lot @Torsten.

[Walter Roberson]

Note: you might want to do something like

N = 5;   %points per interval


tspan = linspace(Delay(i), Delay(i+1), N);

The difference is that when you use exactly two elements in a time range, then ode45() will return back as many time points as it feels like, whereas if you pass in a vector of more than two elements then ode45() will return at those particular times.

[Walter Roberson]

    T = [T ; t(1:end-1)];
    X = [X ; x(1:end-1,:)];

However at the very end, you would want to append t(end) to T and x(end,:) to X .

Also, in some cases you need to adjust the boundary conditions. For example if you were adding a jolt of energy or were bouncing off a floor, then one or more of the boundary conditions for the new interval should reflect the change. "Your y velocity was -5 when you hit the wall, you need to reverse that to +4.97 at the bounce". Because of that, your conditions upon "entering" time Delay(i+1) "from below" might not be the same as your conditions "leaving" that at what is now Delay(i) [because i value changed]. You need to figure out which of the two sets of values you want to put into your accumulated data. @Torsten chose to leave out the "entering" state data and include the "leaving" state data. Which is a perfectly reasonable choice in many circumstances... but sometimes you want the other choice, and sometimes you want the data to include both so that you can visually see the step when you plot.