How to plot dirac n-th derivation ?

I dont know how to plot this function which contains first and second derivation of dirac.

This is how i tried to solve it

Y= 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*dirac(t,1) + 159*dirac(t,2) + 124126*(t^3);

ezplot(Y ,[0 1])

But Matlab send me:

Error using sym/dirac

Too many input arguments

Can someone please help. Thanks

Answer 1

[Star Strider]

You have the arguments to dirac reversed.

From the documentation:

• dirac(n,x) represents the n-th derivative of the Dirac delta function at x.

[joejonson]

Thanks you for the answer.

That is what confusing me, beacuse i get this as answer when finding inverse Laplace of some another function. Do you know maybe how can i plot this ?

[Star Strider]

This works for me:

syms t

Y= 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*dirac(1,t) + 159*dirac(2,t) + 124126*(t^3);

figure(1)
ezplot(Y ,[0 1])

I suggest you send the code that created your original ‘Y’ to MathWorks as a bug report. The Symbolic Math Toolbox should be producing consistent code. It is quite possible that you’re the first to discover this problem. Also, include the URL of this post, http://www.mathworks.com/matlabcentral/answers/238437-how-to-plot-dirac-n-th-derivation, to give it a context.

[joejonson]

It still doesn't work to me, but thanks you very much for helping.

[Walter Roberson]

joejonson, which MATLAB version are you using?

[joejonson]

Matlab R2012b ,why are you asking Walter ?

[Walter Roberson]

The MuPAD version of dirac() has the derivative number after the expression but the interface that the Symbolic Toolkit provides from MATLAB to call MuPAD has the derivative number before the expression. If you attempted to use MuPAD code directly in MATLAB then the arguments would be in the wrong order.

The present order has been used at the MATLAB level since R2007b when MuPAD replaced Maple as the symbolic engine.

Could you check which version of dirac you are invoking?

which -all dirac

[joejonson]

F:\Program Files\MATLAB\R2012b\toolbox\symbolic\symbolic\dirac.m

F:\Program Files\MATLAB\R2012b\toolbox\symbolic\symbolic\@sym\sym.m % sym method

Is that it ?

[Walter Roberson]

That does look appropriate.

If you use Star Strider's code,

syms t
Y = 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*dirac(1,t) + 159*dirac(2,t) + 124126*(t^3);

then do you get the error at that point, or do you get the error when ezplot is processing it? If you use

subs(Y, t, 5.4321)

does that trigger the error? If not then does double() of the result trigger the error ?

[joejonson]

Error occurs beacuse these two dirac (dirac(1,t) + dirac(2,t)) in both case. Subs also trigger error.

[Walter Roberson]

As an experiment, try

syms t
Y = 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*diff(dirac(t,t) + 159*diff(dirac(t),t,t) + 124126*(t^3)

what does it return?

[joejonson]

Error using sym/dirac

Too many input arguments.

Error in dsdsa (line 2)

Y = 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*diff(dirac(t,t)) + 159*diff(dirac(t),t,t) + 124126*(t^3)

[Walter Roberson]

Correction,

syms t
Y = 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*diff(dirac(t),t) + 159*diff(dirac(t),t,t) + 124126*(t^3)

[joejonson]

There is no error now.

This is what I got.

(10078*t)/125 + 731939*dirac(t) + (6373*t^2)/1000 + 124126*t^3 + 4663809*dirac(t, 1)

[Walter Roberson]

I do not know why I didn't find it before, but according to the release notes, the two-input form of dirac is new as of R2014b.

Are you able to ezplot() the version that uses dirac(t,1) ?

It appears that it is treating dirac(t,2) as being equal to dirac(t,1) which does not appear to be correct.

I would be interested to see the result of

evalin(symengine, '80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*diff(dirac(t),t) + 159*diff(dirac(t),t,t) + 124126*(t^3)')

and

evalin(symengine, '80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*dirac(t,1) + 159*dirac(t,2) + 124126*(t^3)')

[joejonson]

Answer is the same:

80.624*t + 731939*dirac(t) + 6.373*t^2 + 124126*t^3 + 4663650*dirac(t, 1) + 159*dirac(t, 2)

[Walter Roberson]

Okay, so the workaround is this:

dirac = @(varargin) feval(symengine, 'dirac', varargin{:})

However, I have a suspicion you might be encountering https://www.mathworks.com/support/bugreports/details/571941 or https://www.mathworks.com/support/bugreports/details/578535, both of which are marked as fixed as of R2010b. Both of them are about problems with simplifying expressions that contain derivatives of dirac. Perhaps the fix did not really "take". Try simplify() the output of one of the evalin() and see if the expression changes.

[joejonson]

Expression is changed, but dirac functions are still the same.

[Star Strider]

Walter, Joe — I apologise for appearing to abandon this. I didn’t intentionally. Had to take a comatose laptop for repairs, an 85 km round-trip.

[joejonson]

It is OK, you dont have to apologise.

Do you maybe have some ideas how to help me with plotting.

[Star Strider]

I thought I did the plot. I don’t have access to R2012b. (This machine has Win 10 so I can’t install anything earlier that R2015a on it; the comatose computer has R2013a.)

With this code:

syms t

Y= 80.624*t + 731939*dirac(t) + 6.373*(t^2) + 4663650*dirac(1,t) + 159*dirac(2,t) + 124126*(t^3);

figure(1)
ezplot(Y ,[0 1])

This is the plot I get:

[joejonson]

Can you please raise time, from 0 to 100,and show me picture. I am expecting some transient in this response,so i am curious.

[Star Strider]

I don’t see any transients:

[joejonson]

Yeah, too bad. Is it problem if you tried it once more time,raising time up to 10000,just in case.. Thanks

[Star Strider]

My pleasure.

No real changes in the shape of the curve, and no transients. I also looked at ranges from 1E-3 to 1E+7.

[joejonson]

Star Strider and Walter Robertson - thanks you very much for everything, both of you really helped me out.

I hope i didnt bother much.

Thanks again!

[Star Strider]

My (our) pleasure. Walter — thanks for picking this up and exploring it while I was out running errands. (Life intrudes.)

No bother, but if you want us to see what you did with your dynamical system to get the result you did, describe it and post your code. No promises, but we might be able to help with it. If you expect transients, something went wrong somewhere. (If I remember correctly, the delta function should integrate to a constant.)

[joejonson]

This is how my code looks like:

 syms s t;
 Rb=0.043;
 Ra=1.071;
 La=6.824;
Lb=0.272;
 M=1.36;
 C=0.000398;
 i1=(200*((Rb + s*Lb + 1/(C*s))/((s^3)*(La*Lb*M - M^2)) + (s^2)*(Ra*Lb*M + La*Rb*M) + s*(M*Ra*Rb + La*M/C) + M*Ra/C));
 F=ilaplace(i1);
 pretty (F)
ezplot(F ,[0 1])
ylabel('I [A]')
xlabel('t [s]')

[Star Strider]

How did you get i1?

Is it a circuit analysis problem? If so, what does the circuit look like?

[joejonson]

It is circuit analysis.This is how it looks like:

There is no much parameters,just transformers primary and secundary resistance (Ra,Rb) and inductance(La,Lb) ,mutual inductance (M) and capacity (C).DC voltage is 200 V (step function).

I get i1 analizing this circuit in Laplace domain (two equations with two unknowns)

[joejonson]

This are the first two equations:

From second equation I excrete I2, and then put it in first equation.

[joejonson]

I wrote the code in a simpler format and now everything is OK. I now got this :

This is something that I was expecting, beacuse at previuos pictures current goes to infinetly, so i was thinking that it will on some point start decreasing to zero,but that does not happend. The problem was created when i want to make final equation simpler, but actually I made it much complicated, and somewhere probably made mistake on multiplying.

[Star Strider]

I’m getting similar results (no transients) with my independent analysis. I can’t find the error. It most likely has to do with my not having done anything with transformers in a very long time.

The Laplace transform of the delta function is 1, and the step is 1/s, and since these are mutliplied by the transfer function, they shouldn’t show up as separate entities in the result.

[joejonson]

I think it is OK now. I was expecting some transient only beacuse current was raising to infinitely. That was strange to me, so I was thinking it will maybe at some point stop raising,and start decreasing to zero. These new results are looking much better.

[Star Strider]

What you’re getting now is what I would expect. The capacitor charges and the current eventually goes to zero in the secondary circuit. I got similar results.

[joejonson]

I am very grateful to both of you, Star Strider and Walter, for your time and all you have done to help me.! I appreciate it!

Thanks you again

[Star Strider]

My (our) pleasure!