My function creates a high oscillating unusual plot

I need to plot and find the optimum point of the function below. But it creates a weird plot which I have inserted. First, I thought this is because the function has coefficients less than matlab eps. But setting them to zero did not change the plot. Also, the optimization alogorithms (fmincon,ga,particle swarm,PSO) gives wrong answers. Does anyone have an explaination for this behaviour?

Thank you very much.

f =17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0

Answer 1

[Walter Roberson]

The summary of the below investigation is that what you are observing appears to be a bug or limitation in MATLAB.

It might be this bug: https://www.mathworks.com/support/bugreports/2422167?ref=ts_R2020b_Update_6 but that is supposed to be fixed by R2021a Update 1.

I will refer this over to Mathworks as it looks like a bug I reported earlier might not be fixed after all.

syms t1
f =17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0

f = 

f2 = collect(f, t1)

f2 = 

fplot(f2, [0 1e-15])

df1 = diff(f2, t1)

df1 = 

fplot(df1, [0 1e-15])

fplot(df1, [0 3e-17]); ylim auto

T1 = linspace(0, 10, 1000).';
f2 = subs(f, t1, T1);
vpa(f2(1:10))

ans = 

plot(T1, double(f2))

[Walter Roberson]

I updated my case with Mathworks.

The summary of the workaround I used here is to choose specific values and subs() them in and then double() the results.

The alternative workaround would be to wrap the function like

F = @(T1) double(subs(f, t1, T1))

[Wathsala Karunarathne]

Thank you very much, Walter.

Answer 2

[Christopher Creutzig]

What you are observing is a limitation of computing with floating point numbers. When computing a value of your function, say, for t1=1e-15, the computer needs to subtract several intermediate results that are close to one another, and that will introduce artefacts.

As Walter has shown, you can alleviate that by doing these numerical calculations in the Symbolic Math Toolbox, where you can use higher-precision numbers than the doubles MATLAB is working in. But you should understand that just pushes the problem further out, and at some point you may need to increase the precision used, by calling digits(30), digits(50), or even higher values. Also, using variable precision arithmetic is much slower than using the double-precision numbers your computer has special hardware for.

syms t1
f = 17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0;
f1 = @(x) subs(f,t1,x); % force MATLAB to use vpa for fplot
fplot(f1,[0 1])

You may or may not be able to use the same workaround when calling optimization functions - but it will cost you a lot of time, and it is probably better to find a numerically more stable form of your input. You may even want to think about approximating it. Here is one way you can do that, using a plot to assess graphically how good the approximation is. I am assuming you already decided to focus on the range from 0 to 1, but the concept should easily generalize:

f2 = series(f,t1,1/2,'Order',5);
f3 = series(f,t1,1/2,'Order',10);
fplot(f1,[0,1],'-.','LineWidth',1.4)
hold on
fplot(f2,[0,1])
fplot(f3,[0,1],'k')
legend(["Original","Order 5 approx","Order 10 approx"])

If, looking at this, you decide the order 10 approximation is a close enough fit in the range you are interested in, using that for the optimization should give you results very quickly.