I managed to change the angle with fixed numbers, but I need it to change over time.

I managed to make the angle change to the specified angle (perturb_magnitude) after a defined time (perturb_time). However, I need the "perturb_magnitude" to change based on a matrix "jj = [ds ddd]" where ds = the time interval at which the angle should change (perturb_time), and ddd = the angle to which it should change (perturb_magnitude). The jj matrix is attached.

Finally, I need the angle to update over time, where both time and angle are contained in the matrix "jj = [ds ddd]." ; "jj = [time angle]."

.
.
.
perturb_magnitude = pi/2;  
perturb_time = 30;  % Time in seconds
perturb_index = find(t >= perturb_time, 1);
x(perturb_index:end, :) = x(perturb_index:end, :) + repmat([perturb_magnitude, perturb_magnitude, perturb_magnitude], numel(t) - perturb_index + 1, 1);

[t_after_perturb, x_after_perturb] = ode45(@(t, x) (equation), t(perturb_index:end), x(perturb_index, :));

Answer 1

[Walter Roberson]

Any time you have a value needed for ode calculation change discretely, you need to stop executing the ode*() call at the boundary of the change. If necessary then adjust the boundary conditions. Then call the ode*() function starting from the time you left off.

If the change in parameter happens at particular times, then this can be handled by looping using tspan that end at the scheduled time of the next change. If the change in parameter is not at predictable times (for example it happens when a rotation angle is detected to cross a particular boundary) then you need to use event functions to signal the stop.

It is common for people to use try to use interp1() inside the ode function. Doing that violates the mathematical conditions for the ode*() call to be valid, unless a non-default interpolation method is used that has continuous second derivatives -- for example you can request cubic spline interpolation and that would satisfy the mathematical constraints. But if you do that, remember that the implication of any interpolation method is that the system is "anticipating" the change and moving towards it, such as starting to rotate so that at the given time it will already reach the known angle. (As compared to a rotation that does not start until a mechanism is tripped.)

[Daniel santiago Umaña]

I tried this, but it didn't work, and I'm also encountering this error.

Error using repmat

Replication factors must be a row vector of integers or integer scalars.

Error in controles6 (line 35)

x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index + 1, 1);

jj = [tSpan ddd];
for i = 1:size(c,1)
    perturb_magnitude(i) = ddd(i);  
    perturb_time(i) = tSpan(i);  
    perturb_index(i) = find(t >= perturb_time(i), 1);
    x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index + 1, 1);
    [t_after_perturb, x_after_perturb] = ode45(@(t, x) (equation), t(perturb_index(i):end), x(perturb_index(i), :));
    t_after_perturb = t_after_perturb(t_after_perturb <= max(tSpans));
    x_after_perturb = x_after_perturb(1:length(t_after_perturb), :);
end

[Walter Roberson]

    x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index + 1, 1);

You are calling repmat() specifying two parameters for replication factors -- numel(t) - perturb_index + 1 in the first parameter, and 1 in the second parameter. repmat() can be called with individual replication parameters only in the case that each of the replication parameters are scalars . But perturb_index is a vector, so numel(t) - perturb_index + 1 is a vector .

Perhaps you wanted numel(t) - perturb_index(i) + 1 ?

[Daniel santiago Umaña]

Yes, thank you, but now I'm facing these issues:

Error using odearguments

When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in controles6 (line 38)

[t_after_perturb, x_after_perturb] = ode45(@(t, x) (A - B * K) * x, t(perturb_index(i):end), x(perturb_index(i), :));

[Walter Roberson]

At some point, perturb_index(i) is coming out the same as the length of t -- so some perturb_time value is greater than the second last element of t but does not exceed the last element of t.

[Daniel santiago Umaña]

Any suggestions on how I can fix this part? As for the programming logic to achieve what I need, do you think it's correct?

Thank you for your attention!

[Walter Roberson]

You need to ensure that you have a final t value that is strictly greater than the last perturbation time, so that t(perturb_index(i):end) will always have at least two elements.

Also, you should be probably be using

[t_after_perturb, x_after_perturb] = ode45(@(t, x) (A - B * K) * x, t(perturb_index(i)+[0 1]), x(perturb_index(i), :));

since you should not be asking ode45 to cross the boundary of the next scheduled perturbation.

[Daniel santiago Umaña]

Thank you!

The other errors have been resolved, but this one appeared:

Index exceeds the number of array elements. Index must not exceed 537.

Error in controles6 (line 37)

[t_after_perturb, x_after_perturb] = ode45(@(t, x) (A - B * K) * x, t(perturb_index(i)+[0 1]), x(perturb_index(i), :));"

[Walter Roberson]

You still need to ensure that the final t value is strictly greater than the last perturbation time, so that perturb_index(i) is less than the index of the last t.

If necessary, choose a duration after the last perturbation time, and set

t(perturb_index(end)+1:end) = [];
t(end+1) = t(end) + duration_after_last_perturbation

This will have the advantage of stopping running after the last perturbation has had time to die out.

But if for some reason you need to run to a fixed time even if nothing is happening, then you may instead need to delete one or more perturbations so that you are not trying to run the ode after the hypothetical fixed end time.

[Daniel santiago Umaña]

Thank you!

I initially tried it like this, increasing the t index by 10 times. The program worked, but it didn't work the way I wanted it to.

for i = 1:size(c,1)
    perturb_magnitude(i) = ddd(i);  
    perturb_time(i) = tSpan(i);  
    perturb_index(i) = find((t*10) >= perturb_time(i), 1);
    x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index(i) + 1, 1);
    [t_after_perturb, x_after_perturb] = ode45(@(t, x) (equation), t(perturb_index(i)+[0 1]), x(perturb_index(i), :));
    t_after_perturb = t_after_perturb(t_after_perturb <= max(tSpan));
    x_after_perturb = x_after_perturb(1:length(t_after_perturb), :);
end

I tried it as you suggested, but it gave me this error:

Arrays have incompatible sizes for this operation.

Error in controles6 (line 39)

x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index(i) + 1, 1);

Related documentation

It turned out like this:

for i = 1:size(c,1)
    perturb_magnitude(i) = ddd(i);  
    perturb_time(i) = tSpan(i);  
    perturb_index(i) = find((t) >= perturb_time(i), 1);
    duration_after_last_perturbation = 10;
    t(perturb_index(end)+1:end) = [];
    t(end+1) = t(end) + duration_after_last_perturbation;
    x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index(i) + 1, 1);
    [t_after_perturb, x_after_perturb] = ode45(@(t, x) (equantion), t(perturb_index(i)+[0 1]), x(perturb_index(i), :));
    t_after_perturb = t_after_perturb(t_after_perturb <= max(tSpan));
    x_after_perturb = x_after_perturb(1:length(t_after_perturb), :);
end

[Walter Roberson]

I am not clear as to whether everything is solved now?

[Daniel santiago Umaña]

Hello, it didn't solve the issue. When I added 'perturb_index(i) = find((t*10) >= perturb_time(i), 1);' the program worked, but not in the way I wanted.

When I added:

't(perturb_index(end)+1:end) = [];

t(end+1) = t(end) + duration_after_last_perturbation'

I encountered this error:

Arrays have incompatible sizes for this operation. Error in controles6 (line 39) x(perturb_index(i):end, :) = x(perturb_index(i):end, :) + repmat([perturb_magnitude(i), perturb_magnitude(i), perturb_magnitude(i)], numel(t) - perturb_index(i) + 1, 1); Related documentation