ezplot fails to plot

I am currently working on a very complicated plot obtained from solving system of 25 differential equations.

The resulting plot should look like an exponential decay. I am 100% sure that the plot range is in the right place.

However, only very short range is plotted in ezplot. It suddenly cuts off at a certain 't' value.

It looks like this:

(exp(44076718174000119488089/4398046511104000000000000000000000000 - (44076718174000119488089*t)/4398046511104)*heaviside(t - 1/1000000000000000000000000)*(4835908122725230768047511886758352784202215283460756517106483200*exp((1377397326799331872275*t)/137438953472 - 55095893071973274891/5497558138880000000000000000000000) - 21822804986935542488305495815699038208*exp((44076450240050567389817*t)/4398046511104 - 44076450240050567389817/4398046511104000000000000000000000000) + 2030859531944623809729207891340025*exp((44076052898326019322713*t)/4398046511104 - 44076052898326019322713/4398046511104000000000000000000000000) - 18085634090625306713632872272494592*exp((44076052898326019322841*t)/4398046511104 - 44076052898326019322841/4398046511104000000000000000000000000) + 899841683603599959456792245698560*exp((44076052898326019323737*t)/4398046511104 - 44076052898326019323737/4398046511104000000000000000000000000) + 16859145488020262757774129136402432*exp((44076444021324750850777*t)/4398046511104 - 44076444021324750850777/4398046511104000000000000000000000000) - 814216110518877199209330571739136*exp((44076499305661593367385*t)/4398046511104 - 8815299861132318673477/879609302220800000000000000000000000) - 59996389385486985376868742957341184*exp((44075613477633556029145*t)/4398046511104 - 8815122695526711205829/879609302220800000000000000000000000) - 5585274059419732529680339483843837206736628160716523821019955200*exp((44076448937248382280633*t)/4398046511104 - 44076448937248382280633/4398046511104000000000000000000000000) + 2237444409888768455238992971956224*exp((44075613477633556027225*t)/4398046511104 - 1763024539105342241089/175921860444160000000000000000000000) + 723346804237244208838146923150757888*exp((44076450624009214494457*t)/4398046511104 - 44076450624009214494457/4398046511104000000000000000000000000) + 2339736958591406137822707047806795776*exp((44076450624009214494649*t)/4398046511104 - 44076450624009214494649/4398046511104000000000000000000000000) - 268336364735555131827998618577920000*exp((44076450624009214494809*t)/4398046511104 - 44076450624009214494809/4398046511104000000000000000000000000) + 499178652185299258833824396018974720*exp((44076450624009214494969*t)/4398046511104 - 44076450624009214494969/4398046511104000000000000000000000000) - 1308598970657934425299226132398012991771653075086868480*exp((8292132738039488089*t)/4398046511104 - 8292132738039488089/4398046511104000000000000000000000000) + 749365938003100155014933091149451847841073721731585070375895040*exp((44074735369117811027289*t)/4398046511104 - 44074735369117811027289/4398046511104000000000000000000000000)))/52656145834278593348959013841835216159447547700274555627155488768

I know this is horrendous work for MATLAB (and my PC), but it is very important that I can plot this because it is a very important data for writing a paper.

IS there any way I can forcefully make MATLAB plot all the way through (at least t = 0 to 0.005)?

[Walter Roberson]

I am waiting for you to reply to http://uk.mathworks.com/matlabcentral/answers/254363-efficient-way-of-calculation-matrix-exponential-in-matlab

[Hayao]

for some reason, I cannot post comments on that thread anymore. I really appreciate all the things you've done for me, but I (as well as my advisors) decided we will not use the results.

I apologize that it has come to this, as I know I made you spend long time on calculations.

Answer 1

[Walter Roberson]

Your sub-expressions that are similar to

exp((44076444021324750850777*t)/4398046511104 - 44076444021324750850777/4398046511104000000000000000000000000)

are overflowing when evaluated numerically with t more than about 0.0005. That is leaving the values undefined.

The function is clearly created symbolically. You will need to keep it symbolic and evaluate it symbolically at an array of specific t, then convert the symbolic result to numeric form to print. When you do that the symbolic exponents can cancel enough to give you a usable result.

[Hayao]

Thank you for the reply!

Well the exponential function obtained from the system of 25 inhomogeneous differential equations is too complicated to solve symbolically. I have to enter a specific value in order to obtain useful data in a reasonable span of time.

Even worse, the obtained function have this "heaviside" functions inside that I am trying to get rid of. Heaviside function arises due to use of Dirac's delta function (as a source of impulse input of the system).

Also, I thank you once again about the matrix exponential question I have made before. I have worked on other systems with similar matrix exponentials, but I was able to solve this through diagonalization of the matrix, which made the calculations significantly simpler and faster.

[Walter Roberson]

You may wish to convert the heaviside to piecewise.

The exponential function might overall be too complicated to solve symbolically, but whatever substitutions you did gave you a symbolic expression that can be evaluated symbolically at particular t values.

In Maple I used

 f0 := (exp(44076718174000119488089/4398046511104000000000000000000000000 - (44076718174000119488089*t)/4398046511104)*heaviside(t - 1/1000000000000000000000000)*(4835908122725230768047511886758352784202215283460756517106483200*exp((1377397326799331872275*t)/137438953472 - 55095893071973274891/5497558138880000000000000000000000) - 21822804986935542488305495815699038208*exp((44076450240050567389817*t)/4398046511104 - 44076450240050567389817/4398046511104000000000000000000000000) + 2030859531944623809729207891340025*exp((44076052898326019322713*t)/4398046511104 - 44076052898326019322713/4398046511104000000000000000000000000) - 18085634090625306713632872272494592*exp((44076052898326019322841*t)/4398046511104 - 44076052898326019322841/4398046511104000000000000000000000000) + 899841683603599959456792245698560*exp((44076052898326019323737*t)/4398046511104 - 44076052898326019323737/4398046511104000000000000000000000000) + 16859145488020262757774129136402432*exp((44076444021324750850777*t)/4398046511104 - 44076444021324750850777/4398046511104000000000000000000000000) - 814216110518877199209330571739136*exp((44076499305661593367385*t)/4398046511104 - 8815299861132318673477/879609302220800000000000000000000000) - 59996389385486985376868742957341184*exp((44075613477633556029145*t)/4398046511104 - 8815122695526711205829/879609302220800000000000000000000000) - 5585274059419732529680339483843837206736628160716523821019955200*exp((44076448937248382280633*t)/4398046511104 - 44076448937248382280633/4398046511104000000000000000000000000) + 2237444409888768455238992971956224*exp((44075613477633556027225*t)/4398046511104 - 1763024539105342241089/175921860444160000000000000000000000) + 723346804237244208838146923150757888*exp((44076450624009214494457*t)/4398046511104 - 44076450624009214494457/4398046511104000000000000000000000000) + 2339736958591406137822707047806795776*exp((44076450624009214494649*t)/4398046511104 - 44076450624009214494649/4398046511104000000000000000000000000) - 268336364735555131827998618577920000*exp((44076450624009214494809*t)/4398046511104 - 44076450624009214494809/4398046511104000000000000000000000000) + 499178652185299258833824396018974720*exp((44076450624009214494969*t)/4398046511104 - 44076450624009214494969/4398046511104000000000000000000000000) - 1308598970657934425299226132398012991771653075086868480*exp((8292132738039488089*t)/4398046511104 - 8292132738039488089/4398046511104000000000000000000000000) + 749365938003100155014933091149451847841073721731585070375895040*exp((44074735369117811027289*t)/4398046511104 - 44074735369117811027289/4398046511104000000000000000000000000)))/52656145834278593348959013841835216159447547700274555627155488768:

Maple uses Heaviside instead of heaviside so it needed to be changed:

f := subs(heaviside = Heaviside, f0):

Then create a sequence of (symbolic) pairs of (t, f(t)), spaced 1/50000 apart and use 250(+1) of them. These gives an upper value of 1/200 = 0.005 at 0/50000, 1/50000, 2/50000, ... 250/50000 .

When particular fractions are substituted in for t, you get terms of the form A*exp(-FRACTION1)*exp(FRACTION2) . Rather than evaluating that numerically, the exponents can be combined to give A*exp(FRACTION2-FRACTION1) that does not have the overflow problems. Maple's combine/expln routine does that conversion

S := combine(([seq])(eval([t, f], t = (1/50000)*T), T = 0 .. 250), expln):

Now you have a sequence of pairs (t,f(t)) that can be converted to numeric form without overflow and plotted in Maple

plots[pointplot](S);

[Hayao]

Okay, thank you. I'll see what I can do.