How to change the heaviside function -- and how to use the new function with symbolic objects?

Hi.

I created this function (it is equal to the heaviside function):

function Y = fzine( X )

Y = zeros(size(X));

Y(X > 0) = 1;

eng = symengine;

if strcmp(eng.kind,'maple')

    Y(X == 0) = nan;

else

    Y(X == 0) = 0;

end

Now, why if

syms x

I get:

>> heaviside(x)

ans =

heaviside(x)

BUT

>> fzine(x) ??? Error using ==> sym.sym>notimplemented at 2621 Function 'gt' is not implemented for MuPAD symbolic objects.

Error in ==> sym.sym>sym.gt at 801 notimplemented('gt');

Error in ==> fzine at 7 Y(X > 0) = 1;

The two functions (fzine and heaviside) share the same code...?

[Andrew Newell]

Note to readers of this question: To see how the accepted answer relates to the question, read How to change the heaviside function -- and how to use the new function with symbolic objects?, How to change the heaviside function -- and how to use the new function with symbolic objects?, and How to change the heaviside function -- and how to use the new function with symbolic objects? in that order.

Answer 1

[condor]

Yes...

and using heaviside is very fast, I wrote also this (very fast):

function Y = fzine( X )
%FZINE Summary of this function goes here
%   Restituisce 0 quando X<0 oppure X==0
%   Restituisce 1 quando X>0
Y = heaviside(X);
if Y==0.5
    Y=0;
end

[Walter Roberson]

If your X is a symbolic object, then heaviside is defined to give you back a symbolic number, which you would have to use double() on in order to compare to 0.5

[Andrew Newell]

The answers are out of sequence now, @Walter. This one is in reply to my suggestion that @condor do the calculation numerically.

[Andrew Newell]

@condor, are we done? If so, it would really help sort out this mess if you accepted the most useful answer and voted for any others that you found useful.

Answer 2

[Andrew Newell]

The symbolic toolbox cannot evaluate X > 0 for a symbolic variable X. You can make a small change in your function so it can handle a symbolic variable:

function Y = fzine( X )

Y = zeros(size(X));

Y(double(X) > 0) = 1;

eng = symengine;

if strcmp(eng.kind,'maple')

    Y(X == 0) = nan;
else

    Y(X == 0) = 0;
end

Now if you try

fzine(sym(1))

you get

ans =

     1

(you can also input a double variable). It still won't handle fzine(x); but notice that

heaviside(x)

just returns the name of the function. It isn't trying to actually evaluate heaviside(x).

[condor]

The problem is that fzine is used in a script where I have a symbolic variable...in that script that is very long I have this structure:

syms x

.....

delta=fzine(x)+....;

.....

y=delta+.....;

solve(y)

So I find "x".

What I don't understand is that using heaviside it works, using a function that has the same code (it changes just in x=0 where my fzine gives 0 and not 0.5) it doesn't work, and I don't understand why....

[Walter Roberson]

fzine = @(x) evalin(symengine,subs('piecewise([x > 0, 1], [Otherwise, 0])'))

[condor]

>> syms x

>> %I try heaviside

>> heaviside(x)

ans =

heaviside(x)

>> %now

>> fzot = @(x) evalin(symengine,subs('piecewise([x > 0, 1], [Otherwise, 0])'))

fzot =

@(x)evalin(symengine,subs('piecewise([x > 0, 1], [Otherwise, 0])'))

>> %BUT

>> fzot(x)

??? Error using ==> mupadmex

mupadmex requires one string argument.

Error in ==> mupadengine.mupadengine>mupadengine.evalin at 98

[res,status] = mupadmex(statement);

Error in ==> @(x)evalin(symengine,subs('piecewise([x > 0, 1],

[Otherwise, 0])'))

>> %DOES NOT WORK.......

[Walter Roberson]

What happens if you try

subs('piecewise([x > 0,1],[Otherwise,0])',x,5)

[condor]

DOES NOT WORK:

>> syms x

>> fzot = @(x) evalin(symengine,subs('piecewise([x > 0,1],[Otherwise,0])',x,5))

fzot =

@(x)evalin(symengine,subs('piecewise([x > 0,1],[Otherwise,0])',x,5))

>> fzot(x)

??? Error using ==> mupadmex

mupadmex requires one string argument.

Error in ==> mupadengine.mupadengine>mupadengine.evalin at 98

[res,status] = mupadmex(statement);

Error in ==> @(x)evalin(symengine,subs('piecewise([x >

0,1],[Otherwise,0])',x,5))

[Walter Roberson]

I mean, what happens if you type that at the command line?

[condor]

>> syms x

>> subs('piecewise([x > 0,1],[Otherwise,0])',x,5)

ans =

1

What does this means?

[Walter Roberson]

It means that piecewise is working. I bet class() of the above result would show sym, a symbolic 1.

So then,

fzine = @(x) subs('piecewise([x > 0,1],[Otherwise,0])','x',x);

[condor]

Yes it is a sym object....

>> fzine = @(x) subs('piecewise([x > 0,1],[Otherwise,0])','x',x);

>> fzine(x)

ans =

piecewise([0 < x, 1], [not 0 < x, 0])

[condor]

Is this working?

[Walter Roberson]

Looks like it.

[condor]

Yes Walter....very very good! It seems working.

Just a question..why do I get 2 results:

>> fzine = @(x) subs('piecewise([x > 0,1],[Otherwise,0])','x',x);

>> syms MSP

>> costi=fzine(MSP-10)*2+fzine(MSP-50)*2.5+fzine(MSP-100)*3+15

costi =

piecewise([100 < MSP, 45/2], [not 10 < MSP, 15], [MSP in Dom::Interval(10, [50]), 17], [MSP in Dom::Interval(50, [100]), 39/2])

>> delta=1.3*costi+costi-MSP

delta =

piecewise([100 < MSP, 207/4 - MSP], [not 10 < MSP, 69/2 - MSP], [MSP in Dom::Interval(10, [50]), 391/10 - MSP], [MSP in Dom::Interval(50, [100]), 897/20 - MSP])

>> double(solve(delta))

ans =

39.1000

34.5000

Answer 3

[Andrew Newell]

Here is a numerical version of what you wish to do. First, the function:

function y = heavy0(x)
y = 0*x;
y(x>0) = 1;

Now solve with it:

costi= @(x) heavy0(x-10)*2+heavy0(x-50)*2.5+heavy0(x-100)*3+15;
delta = @(x) 1.3*costi(x)+costi(x)-x;
fzero(delta,0)

ans =

   39.1000

This also works for problems where the solution is the zero point of a heaviside function:

costi = @(x) heavy0(x-10)+x-10;
fzero(costi,1)

ans =

    10

Answer 4

[Paulo Silva]

good question, I'm also clueless about why that happens

Answer 5

[Walter Roberson]

Your fzine is probably not in the private directory of the symbolic functions, so probably you are not getting the correct methods invoked for all items.

Answer 6

[condor]

So what you mean with private directory (where is it?)... where should I put the fzine?

Answer 7

[condor]

ok I did:

>> which heaviside C:\Program Files\MATLAB\R2010b\toolbox\symbolic\symbolic\heaviside.m

so I cut and pasted the fzine in that folder; now:

>> syms x >> fzine(x) ??? Undefined function or method 'fzine' for input arguments of type 'sym'.

IT doesn't work, help please.

[condor]

Matlab doesn't see the function...

[Walter Roberson]

rehash toolbox

Answer 8

[condor]

>> rehash toolbox >> syms x >> heaviside(x)

ans =

heaviside(x)

>> fzine(x) ??? Error using ==> sym.sym>notimplemented at 2621 Function 'gt' is not implemented for MuPAD symbolic objects.

Error in ==> sym.sym>sym.gt at 801 notimplemented('gt');

Error in ==> fzine at 7 Y(X > 0) = 1;

Still doesn't work...

Answer 9

[condor]

Andrew you say:

"The symbolic toolbox cannot evaluate X > 0 for a symbolic variable X."

Why this doesn;t work with heaviside? It doens't evaluate the function but it accept a symbolic object as input so you can solve an equation that contains heaviside(x)...

??

[condor]

An example:

>> syms x

>> delta=heaviside(x)+15-12*x

delta =

heaviside(x) - 12*x + 15

>> solve(delta)

ans =

4/3

>> delta=fzine(x)+15-12*x

??? Error using ==> sym.sym>notimplemented at 2621

Function 'gt' is not implemented for MuPAD symbolic objects.

Error in ==> sym.sym>sym.gt at 801

notimplemented('gt');

Error in ==> fzine at 8

Y(X > 0) = 1;

Answer 10

[condor]

Here where the error gets generated:

          function X = gt(A,B)
        %GT     Symbolic greater-than.
        A = sym(A);
        B = sym(B);
        if isa(A.s,'maplesym')
            X = A.s > B.s;
        else
            notimplemented('gt');
        end
        end

[condor]

I think the point is that heaviside has something to do with maple symbolic engine; infact also sign(x) with x being a simbolic object doesn't work:

>> syms x

>> sign(x)

??? Undefined function or method 'sign' for input arguments of

type 'sym'.

[Andrew Newell]

Right. Heaviside is calling a function in the MuPad engine.

[Walter Roberson]

Odd... your path says 2010b, which is restricted to using MuPad and cannot use Maple. I wonder if the internals of the symbolic toolbox still say Maple in places ??

[condor]

You are right Walter...and I don't know why! But heaviside (I don't know how) is able to bypass that check above...I mean that with heaviside

isa(A.s,'maplesym') returns 1 otherwise I would get the same error!

If I am able to discover how to bypass that I could replicatye it in my function fzine....

Answer 11

[Andrew Newell]

Here is a different way you could approach this problem:

fzine = (x + abs(x))/2;
delta = fzine+15+12*x;
simplify(solve(delta))

ans =

    -5/4

In fact, you could define your function this way:

function Y = fzine( X )
Y = (X + abs(X))/2;
Y(X~=0) = Y/X;

This now behaves almost exactly like the Heaviside function. Even the derivative works:

>> fd = simplify(diff(fzine(x)))

fd =

    -(abs(x) - x*sign(x))/(2*x^2)
>> subs(fd,-2)

ans =

     0
>> subs(fd,2)

ans =

     0

The only difference is for x=0:

>> subs(fd,0)

ans =

     NaN

(the derivative of heaviside gives Inf).

EDIT: Another difference occurs if you try this:

y = fzine(x);
subs(y,0)

ans = 

    NaN

EDIT: I think the problem with the original function is not just the greater than sign. It is the conditional commands, which make it impossible for the function to return an explicit symbolic expression. For more general cases, you could replace a lot of if/then conditions by expressions involving ... (drum roll) heaviside.

[condor]

Dear Andrew,

I am sorry but this does not resolve any problem...

you are not using the fzine function, you are just resolving an equation in x

[condor]

I need a function that gives

0 when x<=0

1 when x>0

and I need to use this function with x being a symbolic object defined as:

syms x

My question is...Is it possibile to do it?

[Andrew Newell]

Mind you, it is not exactly the same as heaviside, so it's risky to use. Why not use heaviside?

[condor]

If x>0 Ex x=4

Y=4 not 1, isn't it?

[Andrew Newell]

Sorry, I got a bit hasty there (I'm trying to keep up with your comments!). I added a line above.

[condor]

no problem

Answer 12

[condor]

I think the best solution up to know is the one that Walter suggested (it seems to be a little bit slow but that is ok for my use...):

fzine = @(x) subs('piecewise([x > 0,1],[Otherwise,0])','x',x);

however I still don't understand:

1) Why using the fzine defined as above and solving and equation as:

>> syms MSP
>> costi=fzine(MSP-10)*2+fzine(MSP-50)*2.5+fzine(MSP-100)*3+15

costi =

piecewise([100 < MSP, 45/2], [not 10 < MSP, 15], [MSP in Dom::Interval(10, [50]), 17], [MSP in Dom::Interval(50, [100]), 39/2])

>> delta=1.3*costi+costi-MSP

delta =

piecewise([100 < MSP, 207/4 - MSP], [not 10 < MSP, 69/2 - MSP], [MSP in Dom::Interval(10, [50]), 391/10 - MSP], [MSP in Dom::Interval(50, [100]), 897/20 - MSP])

>> double(solve(delta))

ans =

39.1000
34.5000

returns 2 solutions?? The first is the correct one...the second is what we would have if in costi all the fzine are = zero so we would have: costi = 15 and costi*1.3+costi=34.5 BUT MSP cannot be 34.5 otherwise costi would not be 15.... !?!

2) I still don't understand why heaviside works with symbolic objects while another function that uses the same code does not...

[Walter Roberson]

1) I don't know. I used the piecewise version in Maple and it gave only the first (391/10) solution.

2) I do not have the Symbolic Toolbox so I cannot check or test the heaviside source code myself.

Answer 13

[Andrew Newell]

@Condor, I think you are mistaken that the function heaviside uses the same code as the other functions we have considered. heaviside is an interface to MuPAD code, which creates a MuPAD object that can be manipulated according to the rules in MuPAD. For example, operators like diff and int can act on heaviside, while they cannot on Walter's version of fzine. My function behaves a little more like a symbolic object, but I have described some shortcomings.

You can easily create a function that gives the correct numerical answer for a numerical input. But you want it to do more than that. You want it to create that MuPAD object even though we don't even know the properties of MuPAD objects.

It's an interesting exercise seeing where attempts to recreate heaviside fall short, but do we really need to go further with this? Why not just use heaviside?

[condor]

Yes I missed the changed you do on your function above, now I'll test it.

I cannot use heaviside for 2 reasons:

1) I am curious to know how to recreate a function that works as heaviside but gets different value

2) heaviside returns 0.5 for x=0 and this is not good for me

[condor]

Excellent Andrew! I think your solution is the best for what I have to do...it is easy!

However at this point I still wait for somebody that explains me the two points above about heaviside and the two solutions

[condor]

The only thing I am experiencing is that it's very very slow to resolve equation...but at this point I am happy with this solution

[Andrew Newell]

Yes, it is much slower than using heaviside.

[Andrew Newell]

The difference in value at x=0 only matters if your solution turns out to be 10, 50, or 100 in the example you've given. But you could check that these aren't solutions by substituting, e.g., subs(costi,10).

[condor]

I cannot because the function that I am looking for will be used in a GUI with a long long code

[Andrew Newell]

Then you would be better off treating this as a numerical problem in Matlab itself instead of MuPAD. Implementing a numerical version of heaviside is trivial, and you can use FZERO to find the solution.

Answer 14

[Konstantinos]

The way I found so you can do your job is for example:

I wanted to find the Z transform of the following:

syms n z
x2 = n*heaviside(n) + (6-2*n)*heaviside(n-6) + (n-6)*heaviside(n-6);
ztrans(x2, n, z)

and it would give me a wrong value because of the heaviside(n-6). So I fixed it by instead using:

 x2 = n*heaviside(n) + (6-2*n)*heaviside(n-5.5) + (n-6)*heaviside(n-5.5);

Basically I guess using some double inbetween your discrete values will fix all the problems. I could have also changed heaviside(n) to heaviside(n+0.5) but there wasn't a point in this particular example. Hope this helps.