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

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

Answer 1

[Torsten]

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

[Nuno Nunes]

I can't seem to place the EventsFunction anywhere. I mean, if I put it before the menu, I get:

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.

If I put it inside the option selected, I get:

Error: File: MainSimulation.m Line: 30 Column: 13
Function definition are not supported in this context. Functions can only be created as local or nested
functions in code files.

And no sense in putting it in the end.

I really don't know where I should insert this function.

Help?

[Torsten]

Save the event function in your working directoy as "EventsFunction.m" and place it nowhere else.

[Nuno Nunes]

I have done so, but MATLAB said:

Local function name must be different from
the script name.

So I actually named it something else (EventsFunct.m)

However, now I get the

SWITCH expression must be a scalar or a
character vector.

Which, from what I uderstand is due to improper type definition.

I've ran the debugger

dbstop if error

And I've found out that the switch function refers to a 'EventFcn' which is a function_handle with value = @EventsFunct. (switch function is on 'odeevents.m')

But when I'm at the EventsFunct.m file, I don't see what it might be.

My EventsFunct.m file is just this:

% Events Function to allow ODE45 to stop integrating, when 'y^2+v^2 = 0'

run('data.m')

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

Help?

[Torsten]

Yes, remove the command run('data.m') or do it in the Events Function.

And rename the script, not the events function file. The function file should have the same name as the function which is defined in it.

[Nuno Nunes]

I'm terribly sorry, I know I must be coming off as a total noob, but I've just got the chance to work with Matlab, for the first time, a couple of months ago.

On the MainSimulation.m I had that, after the ODE45 ran, it should show me the plot of Y(:,1) and Y(:,2).

So, if my ode45 function is likes this:

[T,Y] = ode45 ('StraightRun', t, y0, opt);

What exactly should I input in the EventsFunction as the squared variables?

Because no matter if I put y(#), Y(#) or y0(#), it always show's me one of these errors:

Unrecognized function or variable 'Y'.

or

Not enough input arguments.

Which both relate to the fact that the variable has not been defined yet.

It doesn't really make sense, because if my code is this:

    if option == 1                                                       % Choosing Straight Run
            Vs = 10 * (1852 / 3600);                                     % [m/s] Ship approach speed (reasonable value)
            a_1 = 309;
            t = (0:a_1);                                                 % Manouevre time
            y0 = [Vs 0 0 0 0 0 0];                                       % Initial condition - surge valocity is equal to the approach speed
            run('data.m');
            run('EventsFunction.m');                                     % This Events Function allows for ODE45 to stop integrating when 'u^2+v^2 = 0'
            opt = odeset('Events',@EventsFunction);                      % ODE Event Location to stop ODE45 when ship stops (Vs = 0)
            [T,Y] = ode45 ('StraightRun', t, y0, opt);                   % Calls function
            

y0 is defined right before running data.m and the ode45.

I've even tried defining the position of the EventsFunction(t,y) as:

position = (y(t,1)^2 + y(t,2)^2); 

To include the 't', but even then it claims

Not enough input arguments.

[Nuno Nunes]

I've been trying to wrap my head around this, but no deal.

Why does the 'Not enough input arguments.', regarding the Events Function appear when running the MainSimulation, when every variable written on the EventsFunction is defined on the .m file, and also every single veriable is given by Straight Run??

Main Simulation:

            Vs = 10 * (1852 / 3600);                                     % [m/s] Ship approach speed (reasonable value)
            a_1 = 309;
            t = (0:a_1);                                                 % Manouevre time
            y0 = [Vs 0 0 0 0 0 0];                                       % Initial condition - surge valocity is equal to the approach speed
            run('data.m');
            run('EventsFunction.m');                                     % This Events Function allows for ODE45 to stop integrating when 'u^2+v^2 = 0'
            opt = odeset('Events',@EventsFunction);                      % ODE Event Location to stop ODE45 when ship stops (Vs = 0)
            [T,Y] = ode45 (@StraightRun, t, y0, opt);                    % Calls function

Events Function:

function [position,isterminal,direction] = EventsFunction(t,y)
  position = (y(t,1)^2 + y(t,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

Error obtained:

Not enough input arguments.

Error in EventsFunction (line 4)
  position = (y(t,1)^2 + y(t,2)^2);                                  % The value that we want to be zero

Error in run (line 91)
evalin('caller', strcat(script, ';'));

Error in MainSimulation (line 32)
            run('EventsFunction.m');                                     % This Events Function allows for ODE45 to stop integrating when 'u^2+v^2 = 0'
 

Help?

[Torsten]

Remove the line

run('EventsFunction.m');                                     % This Events Function allows for ODE45 to stop integrating when 'u^2+v^2 = 0'

The Event function is only called from ode45. You don't need to run it.

I don't know if the command " run('data.m') " will cause problems since I don't know this file.

Further, replace

position = (y(t,1)^2 + y(t,2)^2);

by

position = y(1)^2 + y(2)^2;

y is an array, not a function that could be called with t as argument.

If the problems remain, you should copy the "ballode" example for using an event function and play with it to get a feeling of what is allowed and what is not allowed.

[Nuno Nunes]

Great! Thanks!

That solved THAT issue, but another arose:

How do I find out the timestep (a_1) the function stops so that, when plotting data on MainSimulation, I have the ability to axis([0 a_1 -50 50])?

[Torsten]

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

[Nuno Nunes]

All I get is this:

>> te

te =

     []

>> ye

ye =

     []

>> ie

ie =

     []

how to get the corresponding values?

[Torsten]

Then ode45 reached the final time of integration or stopped integration earlier, but the event didn't happen.

Test t(end) and plot y(:,1).^2+y(:,2).^2 over t.

I think

position = y(1)^2 + y(2)^2;

is a bad definition of an event test for y1 = y2 = 0 because a sign-change in the variable will never happen (y(1)^2 + y(2)^2 will always remain >=0 and will never become < 0).

Better use

position = y(1)^2 + y(2)^2 - (a small number),

e.g. 1e-4.

[Nuno Nunes]

What does this Warning mean? It keeps bugging me

Warning: Failure at t=9.169108e+00.  Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (2.842171e-14) at time t. 
> In ode45 (line 351)

And what can I do to correct it?

[Torsten]

It means that the integrator is not able to solve your differential equations for times greater than t=9.169108e+00.

The possible reasons can be manifold (singularity in your solution, too strong error tolerances, "false" integrator (nonstiff instead of stiff), programming errors etc.).

I suggest you analyze the solution you get up to t = 9.16 and see whether it is reasonable.

[Nuno Nunes]

Is there a way to create an EventsFunction but instead of interrupting the ode45 when the 'position' is equal to 0, make it act when it's below a certain value? Because I've had

  position = y(1)^2 + y(2)^2 - 1e-4

But the thing is, since ode45 is giving discontinuous values, it will only act if ode45 gives out EXACTY "1e-4"

Which, from my experience, it does not.

Help

[Jan]

@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.

Answer 2

[Walter Roberson]

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.

If StraightRun needs additional data, see http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Nuno Nunes]

Data.m is a file that gives the vessel's characteristics and has around 150 lines of code defining the components of the hydrodynamic forces present, which are read by Straight Run.

Straight Run, also has 100 lines of code, but I'd rather keep it separate, because on the menu there is another Manoeuvre to be selected and that's a little heavier. I'd like them to be kept on separate files.

Straight Run (and Zig Zag) also call for another file, which hasn't really given me any headache (so far): wind.m, which is the wind coefficient acting on the vessel (hull and sail), whose information is given in data.m.

y(1) and y(2) are supposed to be the instant velocities in x and y, respectively, but I haven't been able to make it work so far.

[Walter Roberson]

If you are initializing matrices used for the calculation, then see the link about parameterizing functions.

Answer 3

[Jan]

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