Different solutions from the same equations in NDSolve ?

Questions
,
1

Dear community members, I am solving the following set of differential equations with NDSolve. It has been a surprise to me that when solving the problem repeatedly I obtain slightly different solutions. To show this I have repeatedly solve the problem 5 times and plotted the solutions. Any ideas why this may be happening?

DutyCycle[d_, Ts_, t_] := If[0 <= Mod[t, Ts] <= d Ts, 1, 0];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
      st1/L vc[t] + (E1 - (st1 Vf))/L,
                   vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], 
    iL[0] == iL0, vc[0] == vc0} /. {R -> 10, L -> 0.2 10^-3, 
    C1 -> 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 0, nc Ts}], {i, nruns}];
Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}]

Best regards
Jesus Rico

2

Are your equations chaotic? I get different solutions from FindMinimum for complex problems if I run them repeatedly.

3

I did 8 plots, and they were all the same (Mathematica Version 7)
Regards
H. Dolhaine

4

I see two different plots when I run the code using Mathematica 10.2 running on Windows 8.1 - 64.

5

Hi again!
My equations may be chaotic but not for the parameters I am using in this case. And Yes, I got two different plots with mathematica 10.2 but in which run the solution will change is a random event. If they were only two different solutions, why is that??

Best Regards

Jesus Rico

6

What happens if you substitute the two solutions into the differential equations? You could calculate and plot an error term for each equation to see if either of the solutions is correct.

7

Frank
I now which of the two is the most accurate because I have worked out some semi analytical solutions and I really can continue with what I am doing but it has being a little intreging why a sophisticated algorithm such as NDSolve behaves a little randomly. I guess that this little problem may be indeed very important when solving the same system near chaotic operating points.
Thanks

Jesus Rico

8

Jesus, What operating system and what version of Mathematica are you using? That info may help understand the problem since Hans doesn’t see it. Hans, same question to you.

9

Frank
I am running this on a Mac Pro, My operating system is Yosemite but I do not know how to get the version (sorry I am not good at that) I have mathematica 10.0 but I have obtained the same inconsistency in Wolfram Desktop 10.0.1.
Jesus

10

Hi,

I suspect that this is a numerical issue related to the If function. If you plot the DutyCycle function with your parameters over time you get:

Plot[DutyCycle[Td, Ts, t], {t, -nc Ts, nc Ts}]

ScreenShot2015-08-15at20.24.03.png
ScreenShot2015-08-15at20.24.03.png381×231 4.81 KB

As t changes from negative to positive the vector field changes.

Manipulate[VectorPlot[Evaluate[{-((Rs (1 - st1) + (Rd st1))/L) iL - st1/L vc + (E1 - (st1 Vf))/L /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3, t -> d}, st1/C1 iL - 1/(R C1) vc} /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3, t -> d}], {iL, 0, 16}, {vc, 28, 30}], {d, -0.1 nc Ts, 0.1 nc Ts}]

Vectorfield.gif

If you plot the trajectory in the phase plane (the initial point is marked by a red dot)

DutyCycle[d_, Ts_, t_] := If[0 <= Mod[t, Ts] <= d Ts, 1, 0];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
      st1/L vc[t] + (E1 - (st1 Vf))/L, 
    vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[0] == iL0, 
    vc[0] == vc0} /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 0, nc Ts}], {i, nruns}];
GraphicsRow[{ParametricPlot[{Table[{vc[t], iL[t]} /. sol[i], {i, 
      nruns}]}, {t, 0, nc Ts}, AspectRatio -> 1, ImageSize -> Medium, 
   Epilog -> {Red, PointSize[Large], Point[{vc0, iL0}]}], 
  Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}, 
   PlotPoints -> 40, ImageSize -> Medium]}]

you get

ScreenShot2015-08-15at20.33.08.png
ScreenShot2015-08-15at20.33.08.png755×378 21 KB

If you start anywhere near (but numerically close to) time zero, you only obtain one solution. So

DutyCycle[d_, Ts_, t_] := If[0 <= Mod[t, Ts] <= d Ts, 1, 0];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
      st1/L vc[t] + (E1 - (st1 Vf))/L, 
    vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[10^-10] == iL0, 
    vc[10^-10] == vc0} /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 10^-10, nc Ts}], {i, 
   nruns}];
GraphicsRow[{ParametricPlot[{Table[{vc[t], iL[t]} /. sol[i], {i, 
      nruns}]}, {t, 0, nc Ts}, AspectRatio -> 1], 
  Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}, 
   PlotPoints -> 40]}]

and

DutyCycle[d_, Ts_, t_] := If[0 <= Mod[t, Ts] <= d Ts, 1, 0];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
      st1/L vc[t] + (E1 - (st1 Vf))/L, 
    vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[-10^-10] == iL0, 
    vc[-10^-10] == vc0} /. {R -> 10, L -> 0.2 10^-3, 
    C1 -> 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, -10^-10, nc Ts}], {i, 
   nruns}];
GraphicsRow[{ParametricPlot[{Table[{vc[t], iL[t]} /. sol[i], {i, 
      nruns}]}, {t, 0, nc Ts}, AspectRatio -> 1], 
  Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}, 
   PlotPoints -> 40]}]

both give

ScreenShot2015-08-15at20.36.26.png
ScreenShot2015-08-15at20.36.26.png618×264 9.83 KB

Ok, now we can try to define the DutyCycle function using Piecewise, which is recommended in the context of ODES.

DutyCycle[d_, Ts_, t_] := 
  Piecewise[{{1, t < d Ts}, {0, t < d Ts}, {1, 5/4 Ts > t > Ts}, {0, 
     t > 5/4 Ts}}];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
      st1/L vc[t] + (E1 - (st1 Vf))/L, 
    vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[0] == iL0, 
    vc[0] == vc0} /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 0, nc Ts}], {i, nruns}];
GraphicsRow[{ParametricPlot[{Table[{vc[t], iL[t]} /. sol[i], {i, 
      nruns}]}, {t, 0, nc Ts}, AspectRatio -> 1, ImageSize -> Medium, 
   Epilog -> {Red, PointSize[Large], Point[{vc0, iL0}]}], 
  Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}, 
   PlotPoints -> 40, ImageSize -> Medium]}]

This results in

ScreenShot2015-08-15at20.52.13.png
ScreenShot2015-08-15at20.52.13.png748×389 13.2 KB

So not there is only one solution, which I hope coincides with your favourite solution.

Cheers,

Marco

11

I get two solutions on my MacBookAir, running Mathematica 10.2 under OS X Yosemite Version 10.10.5

12

Marco,
Excellent demonstration, big thanks!. I will use the suggested definition as the basis for my periodic DutyCycle function. Also, the solution you have obtained coincides with my semi analytical one.

Best regards
Jesus Rico

13

The code I have that uses FindMinimum and doesn’t give reproducible results for large problems involves a function with a discontinuous first derivative from a Max function. Perhaps discontinuity is the cause of the problem.

14

Hello Frank,

I use Windows 7 and Mathematica 7.0.0. .

Running it 15 times yieds always the same plot:

DutyCycle[d_, Ts_, t_] := 
 If[0 <= Mod[t, Ts] <= d Ts, 1, 0]; E1 = 16; Vf = 0.8; Ts = 
 100 10^-6; Td = 0.25; Rs = 0.001; Rd = 0.001; iL0 = 3; vc0 = 30; \
nruns = 15; nc = 2; st1 = 
 DutyCycle[Td, Ts, 
  t]; Eqns1 = {iL'[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] - 
     st1/L vc[t] + (E1 - (st1 Vf))/L, 
   vc'[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[0] == iL0, 
   vc[0] == vc0} /. {R -> 10, L -> 0.2 10^-3, C1 -> 0.2 10^-3}; Do[
 sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 0, nc Ts}], {i, 
  nruns}]; Plot[{Table[vc[t] /. sol[i], {i, nruns}]}, {t, 0, nc Ts}]

plot-ndsolve.bmp (327 KB)

15

Interestingly, the two different results from the original formulation alternate.

ndsolve.jpg
ndsolve.jpg583×239 16.8 KB

16

…but not on my system :slight_smile:

17

Hi Frank,

I don’t think that in general the solutions alternate.

ScreenShot2015-08-17at19.17.12.png
ScreenShot2015-08-17at19.17.12.png736×268 26.5 KB

Also, the same behaviour occurs on MMA 9. But all of this can be “fixed”, by using the Piecewise function which is recommended for ODEs. If should be avoided in these cases.

Best wishes,

Marco

18

When you went to the Piecewise function, you also got rid of the Mod function inside. Recently I was using CompilePrint to look at a compiled Piecewise function and an equivalent If statement and they looked the same. Could the Mod be the problem?

19

DutyCycle[d_, Ts_, t_] := Piecewise[{{1, 0 <= Mod[t, Ts] <= d Ts}}];
E1 = 16;
Vf = 0.8;
Ts = 100 10^-6;
Td = 0.25;
Rs = 0.001;
Rd = 0.001;
iL0 = 3;
vc0 = 30;
nruns = 5;
nc = 2;
st1 = DutyCycle[Td, Ts, t];
Eqns1 = {iL’[t] == -((Rs (1 - st1) + (Rd st1))/L) iL[t] -
st1/L vc[t] + (E1 - (st1 Vf))/L,
vc’[t] == st1/C1 iL[t] - 1/(R C1) vc[t], iL[0] == iL0,
vc[0] == vc0} /. {R → 10, L → 0.2 10^-3, C1 → 0.2 10^-3};
Do[sol[i] = NDSolve[Eqns1, {iL[t], vc[t]}, {t, 0, nc Ts}], {i, nruns}];
(Plot[{Table[vc[t]/.sol[i],{i,nruns}]},{t,0,nc Ts}])
Table[Plot[vc[t] /. sol[i], {t, 0, nc Ts}], {i, nruns}]

community_ndsolve_problem1.jpg
community_ndsolve_problem1.jpg583×239 17.1 KB

20

Hi Frank,

you are quite right. Mod also seems to cause problems. As you say:

DutyCycle[d_, Ts_, t_] := Piecewise[{{1, 0 <= Mod[t, Ts] <= d Ts}}];

does not seem to work. And

DutyCycle[d_, Ts_, t_] := Ceiling[-(Mod[t, Ts] - Ts d)]

does not work either. So there is a problem with Mod. But

DutyCycle[d_, Ts_, t_] := HeavisideTheta[Sin[2 Pi/Ts*t + Pi*d] - Sin[Pi*d]];

and

DutyCycle[d_, Ts_, t_] := If[Sin[2 Pi/Ts*t + Pi*d] - Sin[Pi*d] > 0, 1, 0];

do work; in spite of the latter having the “If”. Curiously enough,

DutyCycle[d_, Ts_, t_] := HeavisideTheta[SawtoothWave[t/Ts - d] - (1 - d)];

does not work.
In the following case, however, there is no Mod and things still go wrong:

DutyCycle[d_, Ts_, t_] := If[IntervalMemberQ[IntervalUnion @@ Table[Interval[{k Ts, Ts d + k Ts}], {k, 0, 4}], t], 1, 0];

This last one, I couldn’t even make work with Piecewise, i.e.

DutyCycle[d_, Ts_, t_] := Piecewise[{{1, IntervalMemberQ[IntervalUnion @@ Table[Interval[{k Ts, Ts d + k Ts}], {k, 0, 4}], t]}}];

does not work either.

DutyCycle[d_, Ts_, t_] := Evaluate[Piecewise[Flatten[Table[{{1, (k - 1) Ts < t < (k + Td) Ts}, {0, (k + Td) Ts < t < (k + 1) Ts}}, {k, -2,3} ], 1]]];

does work, but

DutyCycle[d_, Ts_, t_] := Piecewise[Flatten[Table[{{If[Evaluate[0 < Mod[t, Ts] < d Ts], 1, 0], (k - 1) Ts < t < (k + d) Ts}, { If[Evaluate[0 < Mod[t, Ts] < d Ts], 1, 0], (k + d) Ts < t < (k + 1) Ts}}, {k, -2, 3} ], 1]];

does not - making the case against Mod even stronger.

If I have not made a mistake, all of these functions should be identical in the interval of interest, say from 0 to nc Ts, which you can test by running

Plot[{DutyCycle[Td, Ts, t]}, {t, -nc Ts, 2 nc Ts}]

— some work some don’t.

Cheers,

Marco

PS: I had some other examples of writing this function, but unfortunately, the session died on me when I wanted to post.