File Exchange Pick of the Week

Functional Programming Constructs 8

Posted by Jiro Doke,

Jiro's pick this week is Functional Programming Constructs by Tucker McClure.

Anonymous functions are extremely powerful and useful, and I use them all the time whenever I want to create a quick function without writing a MATLAB file. For me, it's one of those features which the more I use, the more I appreciate. Loren has written this post on anonymous functions, and she has discussed the topic in other posts as well.

For those of you who are not familiar with anonymous functions, let me give you a simple example. Tapping into my mechanical engineering background, let's simulate a damped harmonic oscillator by solving the second-order differential equation numerically with ode45.

We solve the ODE by passing in the system of first order differential equations in the form of

A second order differential equation becomes a system of two first order differential equations.

function dY = odeFcn(t, y, m, k, b)
  dY = [y(2); -1/m * (k * y(1) + b * y(2))];
end

Then you can solve the system like this (after defining m, k, b, and T):

[~, Y] = ode45(@(t, y) odeFcn(t, y, m, k, b), T, [0, 1]);

Simple enough. However, instead of creating a separate function file for odeFcn, I can just create an anonymous function. I can do this because this is a simple, one-statement function.

% Define mass, spring stiffness, and damping
m = 5; k = 5; b = 1;

% Create anonymous function for the damped harmonic oscillator
odeFcnAnon = @(t, y) [y(2); -1/m * (k * y(1) + b * y(2))];

Now, we can just pass the function handle (to the anonymous function) odeFcnAnon to ode45:

% Define time vector
T = 0:0.1:25;

% Solve
[~, Y] = ode45(odeFcnAnon, T, [0, 1]);
plot(T, Y);
legend('Position', 'Velocity');

Wasn't that great? I didn't have to write a separate file or a function. I just created the ODE function "on the fly". However, I can't create any function on the fly. An anonymous function can only contain one executable statement, i.e. I cannot have multiple statements or flow control commands like if and for.

Well, that's a bummer! I wanted to make my dampled oscillator behave a little differently by adding an external force at a certain time. To do that, my ODE function would look something like this:

function dY = odeFcn(t, y, m, k, b, t0, F)
  if t > t0  % after time t0, apply force F
    dY = [y(2); -1/m * (k * y(1) + b * y(2) - F)];
  else
    dY = [y(2); -1/m * (k * y(1) + b * y(2))];
  end
end

This certainly has more than one statement. How can I possibly do this with an anonymous function? Is my only choice to write a separate function?? ...

Not anymore! My colleague, Tucker, comes to the rescue! He has created a number of functions for flow control commands, as well as other functions that aide in use inside anonymous functions. For example, he has a function iif, which is the functional form of if.

iif(<if this>,      <then that>, ...
    <else if this>, <then that>, ...
    <else>,         <then that>);

Let's see what our ODE anonymous function would look like with this.

odeFcnAnon = @(t, y, t0, F) iif( ...
   t > t0, @() [y(2); -1/m*(k*y(1)+b*y(2)-F)], ... % if t > t0
   true  , @() [y(2); -1/m*(k*y(1)+b*y(2)  )]);    % else

And, let's solve. Notice that I'm passing in two extra arguments, t0 and F.

t0 = 15;   % Force applied at 15 seconds
F = 5;     % External force
[~, Y] = ode45(@(t, y) odeFcnAnon(t, y, t0, F), T, [0, 1]);
plot(T, Y);
legend('Position', 'Velocity');

That was cool. I was able to create an anonymous function for my function that contained conditional statements.

There are lots more to Tucker's functions, and he has written a detailed document deriving some of those functions step by step. It's well-written for such a complex topic. I highly recommend reading through it.

Let me show a couple of more functions, and I'll leave the rest for you to explore. One of the things that are not possible with anonymous functions is the ability to deal with outputs of functions. For example, inside an anonymous function, I am not able to access the second output of a function, or index into the first column of the output. This is because I can't have multiple statements inside an anonymous function, like this:

@() [var1, var2] = myFunction;var2(:, 1)

For this, Tucker has created output and paren. Let's create an anonymous function that would solve the ODE, extract the position vector, and plot it.

plotODEFcn = @(t0, F) plot(T, ...
   paren(...
         output(...
                @() ode45(@(t,y) odeFcnAnon(t,y,t0,F), T, [0, 1]), ...
                2), ...   % Get the 2nd output from ODE45  -  Y
         ':', 1));        % Get the 1st column of Y        -  Y(:, 1)

A little confusing? Don't worry. I was confused as well at first. If you read Tucker's documentation, you'll get a better understanding.

Now, we can easily create multiple plots with different forcing conditions.

figure; hold all;
plotODEFcn(15, 5);   % Force of 5 at 15 seconds
plotODEFcn(10, 5);   % Force of 5 at 10 seconds
plotODEFcn(20, 10);  % Force of 10 at 20 seconds
plotODEFcn(5, 2);    % Force of 2 at 5 seconds
plotODEFcn(0, 0);    % No external force

xlabel('Time'); ylabel('Position');
legend('t0=15, F=5', 't0=10, F=5', 't0=20, F=10', 't0=5, F=2', 'No force', ...
   'Location', 'NorthWest')

Comments

I found Tucker's entry very educational. I learned quite a few techniques I wasn't aware of. If you would like to have some discussion on this topic, feel free to leave comments here or you can visit this post in Loren's blog where Tucker's functions are also highlighted.


Get the MATLAB code

Published with MATLAB® R2012b

8 CommentsOldest to Newest

The function odeFcn has a discontinuity at t0, but Runge-Kutta integrators with stepsize control require smooth functions. Otherwise 3 different behaviours can appear:

1. ODE45 integrates right over the discontinuity without noticing this. Then the result has a poor accuracy only.

2. ODE45 finds the discontinuity and stops with a warning message “Unable to meet integration tolerances without reducing the step size below the smallest value allowed”.

3. ODE45 reduces the step size to such a tiny value, that the integration takes hours to run and the accumulated rounding error dominates the solution. This is exactly what the step size control should avoid.

Which of these 3 effects you see depends critically on the initial values, such that e.g. adding EPS can have a substantial influence although the solution is stable.

In consequence, integrating discontinuous functions is a bad idea from a scientific point of view. Any functions, which cause such discontinuities should be avoided in a function to be integrated: IF, ABS, RAND, NOW, ROUND, MOD, … Because this problem can influence the result substantially caused by numerical effects, they should be seen as conceptional bugs. Even a professional lab for simulating crash tests for a famous German car manufacturer had overseen this problem.

Therefore I’m not happy to see this example here, especially because I warn users in the forum frequently about this topic.

Some references:
http://saba.kntu.ac.ir/eecd/ecourses/Matlab%20premier/Chapter08.pdf
http://comjnl.oxfordjournals.org/content/13/4/401.full.pdf
http://pcmap.unizar.es/numerico/upload/DisconRep-2002.pdf

The last one explains a method to let a DOPRI45 detect and handle a discontinuity. It would be a nice enhancement for Matlab’s ODE45. It could be expanded to handle parameter changes triggered by event functions without the need to restart the integration.

But without doubt, this example shows how to program a discontinuous function clean and clear. Thanks!

Kind regards, Jan

Thanks for pointing that out, Jan! I apologize and am a little embarrassed for my engineering inaccuracy…

In fact, back when I was in graduate school, I was working on a multi, rigid-body dynamics problem with discontinuous forcing functions, and the way I dealt with it was exactly what you mentioned: break up the simulation into the different phases, and simulate each phase separately, using the final state of a phase as the initial condition for the next phase. In my case, I had mechanical stops, so I also accounted for any energy loss that were introduced in the impact, which affected the initial conditions for the next phase.

Yes, if our ODE45 handled this gracefully, that would be great! Have you submitted an enhancement request for this? :)

Not to disparage Tucker’s submission, most if not all IIF (which are actually normally called IFF) forms can easily be converted into a one-liner that Matlab’s anonymous function can accept. In your specific case, for example:

odeFcnAnon = @(t,y) [y(2); -1/m * (k*y(1)+b*y(2)-F*(t>t0))];

Such one-liner functions are generally faster in run-time than their functional equivalent, because of the reduced function-call overhead, and since apparently the JIT/accelerator processes them better than an if-else construct. They are also fully vectorizable, which an if-else construct is not.

The tradeoff here, as in many other performance hacks, is often reduced maintainability and readability. But in the specific example above I dare say that the one-liner is actually more readable/maintainable than the more verbose version…

@Yair,

Thanks for your solution. That does seem more elegant, and I see your point about performance and vectorization.

Perhaps my example didn’t do justice to Tucker’s IIF. I think the Inline Conditionals section in Tucker’s document shows a more pertinent example that highlights the use case of a functional IF.

Anonymous functions as a concept are ok, but the syntax is horrible. The worst thing mathworks has done is to hire computer science graduates. They have apparently infiltrated the core of matlab and made it depart from the consistency of the language. The founder of matlab was great in his concept and vision – unfortunately in the last 5-7 years I have witnessed horrible intrusions into the consistency of the language. I knew things were going downhill when the ‘code checker’ suggested me for the first time that strvcat is not cool. Anonymous functions, and in particular function handles (their syntax and concept all together) were further signs that mathworks is run by less than competent scientists – the main target audience for matlab. I don’t need bookkeepers to tell me what to do in matlab – I need smart scientific-oriented leaders to shape the future of matlab – the platform so great in its initial concept. Most importantly, please keep it consistent. I don’t need some ‘dont know better’ to throw ‘objects’, ‘methods’ and similar at me. Give me functions, left hand side out, right hand side in. And please don’t pretend as ‘methods’ are in a way divine, or that we ‘dont’ get it’. The best solution is, please remove all computer ‘science’ majors from mathworks, and hire competent scientific people instead.

Yair,

Be careful when using multiplication by logicals as masking. The expression:

y = zeros(size(x), class(x));
y(x>0) = x(x>0);

is equivalent to:

y = x.*(x>0);

except when x contains nonfinite values (where the multiplication could change the value) or is of an integer class (in which case it will simply error; integers cam be multiplied by other integers of the same class or by scalar doubles. Logicals are neither of those.) Consider both sections of code with one of these x vectors:

x = [1:3 -Inf 5 NaN 7:10];
x = int8(-5:5);

In the ODE solver case, having one of the inputs to the ODE function be nonfinite is generally a sign that you’re already in trouble and the ODE solvers don’t support integer data. But in a more general case, these might be perfectly appropriate inputs.

This is a great post, submission and topic! Thank you to everyone involved. When I got back into matlab recently, I was looking into functional programming support and at the time it seemed that matlab didn’t leverage tail recursion to limit the size of the function stack where applicable – is this still the case and, if so, any chance this will change in the future? Cheers!

Jan,

Great comments! I’m always happy when people understand the particulars of numerical integration. I would add that integrating over discontinuities, however, is common in simulation. E.g., suppose something switches on in the simulation, instantly effecting acceleration. The good news is that this can be handled responsibly! The odeset allows a user to specify a time for “events” (discontinuities), causing the integrator to step to t- and t+ (where t is the time of an event), and then continue with its normal time steps. This makes realistic simulation possible. Simulink also handles this beautifully by using zero-crossing detection.

These postings are the author's and don't necessarily represent the opinions of MathWorks.