Where/How exactly do I create an ODE Event Location Function?

16 views (last 30 days)
Hello all,
Here's the thing:
I have a main file (MainSimulation.m), which basically creates a menu in which I can choose a type of marine Manoeuvre to perform, given variables in other two files (data.m and wind.m), which are the vessel-related parameters and the wind conditions + wind coefficient (forces/moment) calculations.
Now, when I choose a Manoeuvre, I must run ode45 to be able to calculate the overall accelerations (F=m.a), like so:
option = menu('Choose a Manoeuvre:','Straight Run','Zig-Zag Manoeuvre','Exit Menu');
%% Straight Run
if option == 1 % Choosing Straight Run
Vs = 10 * (1852 / 3600); % [m/s]-Ship approach speed (reasonable value)
a_1 = 300;
t = (0:a_1); % Manouevre time
y0 = [Vs 0 0 0 0 0 0]; % Initial condition - surge valocity is equal to the approach speed
[T,Y] = ode45 ('StraightRun', t, y0); % Calls function
drift_angle = -atan(Y(:,2) ./ Y(:,1));
In which Y will be a 7 collumn double array [Y(:,1) = speed in X; Y(:,2) = speed in Y; the other collumns are irrelevant rn]
Well the thing is, given the fact that this is a wind-powered vessel (obtained from data.m) I have come to the conclusion that there is a problem with the ODE function (any of them) when the speed (Vs) of the vessel comes to 0, because there are several parameters (regarding the definition of forces/moments) that directly depend on the vessel's speed, and the following Warning appears, if the ODE45 is not done running, given t (which I assume to be due to the vessel coming to a full stop and there still being some integration steps left):
Warning: Failure at t=3.091528e+02. Unable to meet integration tolerances without reducing
the step size below the smallest value allowed (9.094947e-13) at time t.
> In ode45 (line 351)
Therefore, I wish to run the ODE45 function, as stated, until Vs = 0.
I have seen that there is a thing called ODE Event Location which allows for, if a certain variable/function comes to 0, the ODE will stop integrating, and return the values already calculated.
However, I'm not really sure WHERE to insert the following function [and I'm a little confused as to what "isterminal" and "direction" are]:
function [position,isterminal,direction] = EventsFunction(t,y)
position = sqrt(y(1)^2+y(2)^2); % The value that we want to be zero
isterminal = 0; % Halt integration
direction = 0; % The zero can be approached from either direction
end
Should I insert the above function (EventsFunction) in the same file as the menu is, prior to it, obviously, or should I insert it inside the 'StraightRun' file?
Because I have tried inserting the function "EventsFunction" just before calling ode45 (after the definition of y0), but MATLAB gives me the following error:
Error: File: MainSimulation.m Line: 32 Column: 1
Function definition are not supported in this context. Functions can only be created as local
or nested functions in code files.
And if I insert the function BEFORE/outside the option = menu (....), the following error appears:
Error: File: MainSimulation.m Line: 29 Column: 1
Function definitions in a script must appear at the end of the file.
Move all statements after the "EventsFunction" function definition to before the first local
function definition.
ANY Help provided will be MUCH appreciated!
Nuno Nunes

Accepted Answer

Torsten
Torsten on 13 Aug 2022
if option == 1 % Choosing Straight Run
Vs = 10 * (1852 / 3600); % [m/s]-Ship approach speed (reasonable value)
a_1 = 300;
t = (0:a_1); % Manouevre time
y0 = [Vs 0 0 0 0 0 0]; % Initial condition - surge valocity is equal to the approach speed
opt = odeset('Events',@EventsFunction);
[T,Y] = ode45 ('StraightRun', t, y0, opt); % Calls function
drift_angle = -atan(Y(:,2) ./ Y(:,1));
and place the function where MATLAB can access it.
function [position,isterminal,direction] = EventsFunction(t,y)
position = y(1)^2+y(2)^2; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
  15 Comments
Jan
Jan on 31 Aug 2022
@Nuno Nunes: What does "ode45 is giving discontinuous values" mean?
If the value of an event function is below 0 in one step and above 0 in the next step, ODE45 interpolates the step to find the time, where 0 is hit exactly.
"from my experience, it does not." - please post some code to explain this. Of course the show event function is valid and should be detected reliably. Maybe reducing the maximum step size helps, if the event function is > 0 at both steps, but was < 0 between the time points.

Sign in to comment.

More Answers (2)

Walter Roberson
Walter Roberson on 14 Aug 2022
You are using
[T,Y] = ode45 ('StraightRun', t, y0, opt); % Calls function
specifying the function to be executed as a character vector. We recommend you stop doing that and start using function handles. Besides being more efficient, function handles are able to refer to local functions defined in the same file.
So you can create a single file
option = menu('Choose a Manoeuvre:','Straight Run','Zig-Zag Manoeuvre','Exit Menu');
%% Straight Run
if option == 1 % Choosing Straight Run
Vs = 10 * (1852 / 3600); % [m/s]-Ship approach speed (reasonable value)
a_1 = 300;
t = 0:a_1; % Manouevre time
y0 = [Vs 0 0 0 0 0 0]; % Initial condition - surge valocity is equal to the approach speed
opt = odeset('Events',@EventsFunction);
[T,Y] = ode45 (@StraightRun, t, y0, opt); % Calls function
drift_angle = -atan(Y(:,2) ./ Y(:,1));
end
function [position,isterminal,direction] = EventsFunction(t,y)
position = y(1)^2+y(2)^2; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
function dy = StraightRun(t, y)
whatever appropriate
end
I have no idea what your data.m is about. You are not defining anything that needs to be passed into StraightRun.
  2 Comments
Walter Roberson
Walter Roberson on 15 Aug 2022
If you are initializing matrices used for the calculation, then see the link about parameterizing functions.

Sign in to comment.


Jan
Jan on 15 Aug 2022
As Walter has written already, use:
[T,Y] = ode45 (@StraightRun, t, y0);
% Instead of the old style:
% [T,Y] = ode45 ('StraightRun', t, y0);
Providing the function as CHAR vector is deprecated for 20 years now and "modern" function handles are the better choice.
Then your problem does not concern only the event function, but all functions. The rules are:
  • A function can be stored in a specific m-file, which ist stored in a folder inside Matlab's path. See doc addpath and doc pathtool.
  • A function can by stored as local subfunction inside an m file. If this file is a "script", the functions must be appended at the bottom. Scripts are m-files, which do not start with the keyword "function".
  • Nested functions are stored inside another function. They share variables with the enclosing function.
Examples:
% m-file: myScript.m
disp('This is a script');
a = myFunc(2);
% Local subfunction:
function out = myFunc(in)
out = in * 3;
end
% m-file: myFunction.m
function myFunction % A function without inputs
disp('This is a script');
a = myFunc(2);
end % <- required to close the function, if any other function uses "end" also
% Local subfunction - only visible inside this file:
function out = myFunc(in)
out = in * 3;
end
% m-file: myFunction.m
function myFunction % A function without inputs
disp('This is a script');
a = myFunc(2);
% Nested function - only visible in enclosing function
function out = myFunc(in)
out = in * 3;
end
end % <- required to close the function, if any other function uses "end" also
% m-file: myFunction.m
function myFunction % A function without inputs
disp('This is a script');
a = myFunc(2);
end
% Another m-file: myFunc.m - visible to all other functions and scripts
function out = myFunc(in)
out = in * 3;
end

Community Treasure Hunt

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

Start Hunting!

Translated by