Piecewise functions in matlab

Not the most difficult thing to do by any means. But here’s a handy conversion from a math formula to matlab. Say you have the piecewise polynomial, m, defined as:


        / 0              if x < 0,
m(x) = |  -2x²(x - 3/2)  if 0 ≤ x < 1,
       |  1 + (x - 1)    if 1 ≤ x < 3/2,
        \ x - 1/4        if x ≥ 3/2.

So in matlab if I have some variable x with some values, i.e.

x = -1:0.01:2;

then I can formulate the above as:


m = ...
  (0               ) .* (x < 0           ) + ...
  (-2*x.^2.*(x-3/2)) .* (0 <= x & x < 1  ) + ...
  (1+(x-1).^2      ) .* (1 <= x & x < 3/2) + ...
  (x-1/4           ) .* (x >= 3/2        );

Now, I can show a plot:


plot(x,m);

piecewise polynomial plot matlab

This works because the conditions in matlab are now logicals that return a vector the same size as x, with 1′s if the condition was true and 0′s otherwise. Multiplied against the value for the condition and added to the next gives the correct solution. I guess this method is somewhat risky in the sense that if you mess up your logicals or inequalities the addition could sum up erroneous values without recognizing the error. But I think the presentation is very nice and is easily broken up if need be.

Note: The polynomial above is from Higher Order Barycentric Coordinates by Torsten Langer and Hans-Peter Seidel.

Tags: , , ,

41 Responses to “Piecewise functions in matlab”

  1. Mollier says:

    Very nice example, thank you for sharing.

  2. m says:

    you are a fucking genius. Thank you so much. I think this will work.

  3. doug says:

    Alec, would you happen to know how to sum an iterated piecewise function? Say f(x) = x+x^2 and you need to sum on an interval [a,b] the iteration of f as long as x_i is in [a,b] else the sum is zero. Thank you in adavance.

  4. ajx says:

    I’m not sure what you mean. Do you mean something like this:

    
    a = 0;
    b = 1;
    f = a:0.1:b;
    max_iterations = 10;
    for i = 1:max_iterations
    f = (f + f.^2) .* (f >= a & f <=b) + ...
           (0) .* (f<a | f >b);
    end
    

    Otherwise what do you mean by “sum an iterated piecewise function” can you explain how you would write this function with math notation?

  5. Tom says:

    it doesn’t work :(
    but i used if/else conditions instead

  6. ajx says:

    @Tom Which doesn’t work?!?!

  7. mac says:

    This didn’t work for my piecewise function :(

  8. ajx says:

    @mac What function are you trying to do?

  9. Joe says:

    THANK YOU SO MUCH!!!!

  10. Asma says:

    Many thanks , it’s very simple and useful :)

  11. john says:

    using your same function m(x) say now i want to plot m(-x+1) how would i do so?

  12. ajx says:

    m’s only a function of x, right before it’s evaluated. So to do what you want you could try:

    
    y = -x+1;
    m = ...
      (0               ) .* (y < 0           ) + ...
      (-2*y.^2.*(y-3/2)) .* (0 <= y & y < 1  ) + ...
      (1+(y-1).^2      ) .* (1 <= y & y < 3/2) + ...
      (y-1/4           ) .* (y >= 3/2        );
    plot(x,m);
    
  13. Joao Borges says:

    Just what I was looking for…
    Many thanks.

  14. Raul says:

    Awfully useful, thank you very much.

  15. John says:

    sir, thank you you are a genius

  16. andre says:

    Hi!
    Please, how to calculate de integral of this function on interval [-1,2] ?

  17. ajx says:

    Andre, remember that you can just calculate the sum of the integrals of each subfunction on its respective interval. A more detailed explanation.

  18. Andre says:

    Dear Ajx,
    Thank you for your answer, but I really would like to know if there is a way to make it directly, that is, once we have the expression, is there a command that make the calculation whitout I have to command interval by interval the respective expression?
    Sorry for my english … my corrent language is portuguese.

  19. ajx says:

    I’m not sure exactly what you want to do. You could try MATLAB’s symlink, or maple or mathematic to compute the integral automatically. All of them should be able to easily handle piecewise functions.

  20. Muhammad Hammad says:

    thanks alot deaR….exactly what I was googling for :)

  21. Frank says:

    How would you do a multi-variable piecewise function?

    F(x,y)=
    0 when x<0 y<0
    x*y 0<=x<=1 0<=y<=1
    x 0<=x1
    y 0<=y1
    1 x>1 y>1

    I tried the method given but it adds up the values
    z= (0) .*(x<0 | y<0)+…
    (x.*y) .*(0<=x<=1 | 0<=y<=1)+…
    (x) .*(0<=x1)+…
    (y) .*(0<=y1)+…
    (1) .*(x>1 | y>1);

    any help would be great, thanks

  22. ajx says:

    Hi Frank,
    I’m not sure I understand. How many variables are there? Just x and y? What are x1 and y1? Also does “x<0 y<0″ mean x less than zero OR y less than zero or does it mean x less than zero AND y less than zero? I have a feeling the component-wise ‘or’ operator ‘|’ is the culprit.

  23. Arcly says:

    Sweeeeet

  24. Frank says:

    / 0 if x<0 and y<0
    | x*y if 0<=x<=1 and 0<=y<=1
    F(x,y)= | x if 0<=x1
    | y if 0<=y1
    \ 1 if x>1 and y>1

    Sorry about that, I must have typed my previous comment half asleep when I was frustrated trying to figure out how to use your method to plot the function. The one I wrote is the correct function, I also changed the “|” to “&”. eventually I plotted the function, though it required a bit of loops. it would be great if there was a way to do it your way. There are 2 variables, x and y.

  25. Abraham says:

    Hi ajx,

    I highly appreciate the fact that you took your time to write such a nice program for graphing piecewise functions. I have tried your code to some other similar functions like the one you have, and it works very well, except for the one that I have to do regarding my project.
    I am working on an estimation problem in wish I need to generate a plot of the following function ( a cubic spline K0 = -3):

    B0(x) =
    0.25*y.^3 if k0 ≤ x < k1,

    -0.75*y.^3+0.75*X.^2+0.75*y+0.75 if k1 ≤ x < k2,

    0.75*y.^3 – 0.75*y.^2 +1 if k2 ≤ x < k3,

    -0.25*y.^3+0.75*y.^2-0.75*y+0.75 if k3 ≤ x < k4,

    where x is generated from
    N =100;
    x = rand(1,N); and decomposed into an interval kx given:

    kx = floor(K*x) and
    y = K*x-kx

    Since B0 has 4 intervals, K = 4.

    The peremeter to be estimated is a sum given by:

    g = aK0*BK0(x) +aK0+1*BK0+1(x) +…aK-1BK-1(x)

    where the other B-splines are translations of B0(x) given by:

    BJ(x) = B0(x-kj). for j = K0, K0+1,….K-1, I also would like to plot these translated version of B0(x). K0 = -3 for cubic splines, and is not the same as for the intervals for B0, k0.

    I would really be glad if you can help me.

  26. MaxPayne says:

    Hi,
    does anybody know how to handle relational operation in symbolic computation. For eg. if I want to check certain condition, if a>=b else NaN.
    does anybody have idea ingerneal about looping in symobilc computation..
    thanks.

  27. kumar sharma says:

    This was a Great Trick. hatsoff :)

  28. Gustav says:

    Awesome!
    Could have read documentation. But thanks to google and you I got to the stuff i needed without reading boring definitions of this and that!

  29. Jake says:

    This is a nice trick, but it has a serious flaw that comes up regularly when I try to use. Consider this example:

    f=@(x)(x>0).*(1./x)+(x<=0).*x;

    When x==0, the function "should" return 0, but it returns NaN. Is there a real way to deal with piecewise functions in matlab?

  30. ajx says:

    @Jake,
    I’ve run into this too. Usually I give up on having a one-liner in these situations, but it seems you can taylor-make something for your situation. There are a few options here: http://www.mathworks.com/matlabcentral/newsreader/view_thread/147044. I like this one for example:

    
    % Assumes x is in ascending order
    g = @(x) [x(x<=0) 1./x(x>0)];
    
  31. Matt says:

    Absolute genius, worked like a charm. Thanks for your help!

  32. Beth says:

    THANK YOU!!!! I’ve been trying to figure out how to do this!!

  33. frank says:

    Hey just curious, I tried using trig such as
    ff = …
    (sin(x)) .* (0 <= x <= 10) + …
    (2) .* (10 <= x);
    plot(x,ff)

    When I use trig it seems to behave differently, any help?

  34. frank says:

    Actually you can delete it haha, I see my mistake. Thanks for the code!

  35. Preet says:

    I want to use truncated ramp function defined as below:
    f(t)=t for t=[0:1:100];
    f(t)=1 for t=[100:1:6000];

    I want to use this function in the form such that I can use f(t+3) + f(t-2) etc.
    How should I define this function?

  36. ajx says:

    @Preet, aside from the fact that your f is defined twice for 100 (does f(100) equal 100 or 1?) then you would use f = @(t) (t<100).*t + (t>=100).*1;

  37. I know this is out of date, but the way I deal with singularities is to make a second variable then substitute that in. It is two lines, but it allows anonymous functions to still be used->

    NewX = @(x) (x==0).*(1) + (x~=0).*(x) 
    MySinc = @(x) (x==0).*(1) + (x~=0).*(sin(pi*x)./(pi*NewX(x)))
    
  38. Not sure where my dot-times went, tho :)

  39. ajx says:

    @Michael, nice. The asterisks got lost during the markdown parsing. I fixed your comment accordingly.

  40. Sanaz says:

    Thanks so musch
    it worked!
    best

  41. Use these proxies to encrypt your entire internet connection thereby making yourself secret.

Leave a Reply