Odesolve bug? Another bug?

[Cornel]

Hi,

Change variable a from below and look at the V2(t)...Why odesolve block is behaving like this? Is it normal?
Prime 10 file attached.
The expected behaviour is like a =7.

[Cornel]

The above results are for this:

IF I change to Fixed:

If I change to Adaptive:

Endless running:

And If I change to Radau:

From a=4 until a=8 its ok, but when a=9 not:

a = 10, 11 - ok

a = 12 - NOT ok...

It seems that Radau doesn't make so many misses (at least in this example), and at least as not many as Adams/BDF does. 

Its a bug? Or what's happening?

[AlfredFlaßhaar]

The solution to this relatively simple initial value problem fails due to the discontinuity of the inhomogeneity. This means that "theoretically" there is no statement about the existence and uniqueness of the solution. Perhaps reformulating it as an integral equation will help. Numerically, less sensitivity can then be expected. Or the inhomogeneity can be approximated by a continuous function. However, MC would probably be unable to cope with this.

[Cornel]

But we do not want to stay with differential eq form (derivative), and not to go to intergrals where odesolve block cannot be used?

Maybe Odesolve is not so good for this problem? I am not sure how could this problem can be translated to other methods of solving differential equations...maybe AdamsBDF Mathcad algorithm? But I do not know at this moment how to translate the above diff eq to others algorithm like AdamsBDF or others...

[ttokoro]

[Cornel]

We should maintain the values for C1 and R1 as there were defined initially, not to changed them...as we like..

[Cornel]

But anyway, the problem is still present even so:
Below is with Adams/BDF:

And with Radau:

For all a:0...12 the same results are maintained in this example with Radau:

But as I said previously, we should not changed the values of components. We should have the initial values:

[Cornel]

Hm, but I see that having Radau option selected:

And with initial condition set to 0V, not 10^(-a), then the things seems to be good so far for this example:

So, the Radau is better than Adams/BDF?

[Cornel]

I found a case where Radau misses:

[Cornel]

With Radau:

With Adams/BDF:

[AlfredFlaßhaar]

As already mentioned, this task cannot in principle be solved over the entire time axis in one go. But it would be possible to proceed step by step according to the levels/intervals of the inhomogeneity (where it is constant in sections). To do this, the solution is calculated up to the end of the first section using the original initial value and then the next constant section is calculated using the final value as the new initial value. But my programming skills are not sufficient for that. In any case, the curves calculated so far over the entire time length with only one initial value should be treated with caution, as there may be peculiarities of the numerics involved.

[Cornel]

What is the first section and what is the second section? Make a drawing.

[Cornel]

So, the question is this: the fault is from Mathcad ode algorithm or the fault is not the fault and in the fact this should be the expected behaviour to misses some pulses some time to time?

[AlfredFlaßhaar]

Sections are the intervals of length, e.g. tend/1000 on the time axis, in which V1 and I(Vx) are constant. The first section starts at t=0 and ends at (as far as can be seen in the diagram) up to t=0.605. The second section starts at t=0.605 and goes up to t=0.61, etc. The solution in the first section has the initial value given at t=0. For the solution in the second section, the ordinate in the final value from the first section is taken as the initial value at t=0.605, etc. The overall solution is therefore put together section by section. This should be clear even without a drawing ;-). I cannot see any other numerical MC-compatible way.

[Cornel]

@Werner_E  what are your thoughts about what you see in this thread? Do you think that this reason described above by Alfred to be the reason on why the V2(t) signal misses pulses?

[Cornel]

Or do You think is rather because Mathcad ode algoritm fail?

[Werner_E]

I have no idea. Maybe you are pushing Prime to its numerical limits.

[Cornel]

But maybe:
@Werner_E could you try this exercise in Mathcad 15 to see if Mathcad 15 deliver same result as Prime 10?
@LucMeekes  could you try this exercise in Mathcad 11 to see if Mathcad 11 deliver same result as Prime 10?

[Werner_E]

Mathcad 15 won't allow units with odesolve

I would suggest defining V1(t) that way

There is an inconsistency in your function I(Vx). its not continuous at Vx=0.7V
The signal jumps from 0.207 uA down to 0A

[Cornel]

Even though, this behaviour is still maintained. Try without units in Mathcad 15. For me my question is if this behaviour is the expected behaviour or not (meaning that there is a fault in Mathcad odesolve algorithm)

This result is with:

This result is with:

But if we change R1 from 10 to 12 then also Radau misses pulses:

[Werner_E]

I did not expect that the suggested definition of V1 would change the behaviour of the sheet. I just suggested it because its much shorter 😉

But in my opinion, however, the discontinuity in I(Vx) should be definitely fixed!

I also don't think that the algorithm(s) used for odesolve is/are faulty but rather it seems to be a 'natural' limitation of a numerical approach, implemented using numbers stored in standard IEEE format.

I am not sure if I could be bothered to retype anything in Mathcad just to find out if we experience similar effects in Mathcad, too, or not. There would be no benefit knowing that and my guess is that we would also experience such effects in Mathcad, too.

[AlfredFlaßhaar]

Please excuse me if I'm butting in again:
The numerical MC solution algorithms are very likely to be error-free if the theoretical background can be assumed to be verifiable. In the present case, this is not the case. The explicit ODE to be solved is V2'=f(t,V1, V2). However, according to the Picard-Lindelöf theorem, the function f must be continuous and a solution only exists in the continuity region. This is not the case in the present case because of the pointwise discontinuity of V! und I(V1). The MC algorithm then calculates something uncontrollably contestable. I see only the section-by-section process described above as a possibility for a numerical solution method in MC. Since the inhomogeneities are section-by-section constant, a shifted exp arc in each section should lead to a sawtooth curve overall due to the negative real eigenvalue of the characteristic polynomial. Unfortunately, my programming skills are not sufficient for this.

[Werner_E]

As @ttokoro  had shown, the numeric solution seems to work pretty well using different input values (without the sheet of course we can only guess that he actually did)

@ttokoro uses a capacity 10000 times larger than originally uses by @Cornel  and so of course a 10000 times larger time constant tau of 5 microseconds (and not 500 ps as in Cornels original approach).

See Re: Odesolve bug? Another bug?
Unfortunately @ttokoro did not post his sheet (as usual) and @Cornel had troubles duplicating what @ttokoro's picture shows. Maybe @ttokoro  did not use a true periodic signal, we don't see which value he had chosen for "a", ...- we don't know.

But do you think that the result (using the larger capacity value) would look significally differently when calculated section-by-section as you suggest from what can be seen in @ttokoro s picture?

[AlfredFlaßhaar]

I have already written something about ttokoro's "solution". In my opinion, there is nothing to add to that. And since his solution was not published, it is not a solution.

Cornel's numerical solution attempt does not fail because of the parameters. Rather, it is fundamentally wrong because the prerequisites (as described, continuity) for the existence of the solution with the impulse approach over the full time are not met. But this is required for the calculation algorithm, and it is based on it.

[Werner_E]

The numerical approximation algorithms are not concerned IMHO as you may think of the square wave signal of being a trapezoid signal , with very steep slopes. Now V1 can be considered continuous and we can expect a solution.
@Cornel  may try modelling his V1(t) as such a trapezoid signal and see if it helps in any way.

I am more concerned about the function I(V2(t)) which (apart from the already mentioned discontinuity which seems to be a modelling error of Cornel) somehow feeds back to the function V2 we are looking for and limits/controls the current in the circuit.

This function I() makes it quite difficult to solve symbolically even for a small range of t where V1 is constant throughout.

[Cornel]

---

@Werner_E wrote:

The numerical approximation algorithms are not concerned IMHO as you may think of the square wave signal of being a trapezoid signal , with very steep slopes. Now V1 can be considered continuous and we can expect a solution.
@Cornel  may try modelling his V1(t) as such a trapezoid signal and see if it helps in any way.

---

 I have made the calculations with V1(t) more as a trapezoid signal:
With:

And with Radau:

[Cornel]

@Werner_E @LucMeekes @AlfredFlaßhaar 
Hm...maybe its has something to do with this, look. So? A true pulse/square wave signal causes problems?

But...which is good:

[Cornel]

[AlfredFlaßhaar]

It is obvious to make the vertical edges/flanks of the individual pulses slightly slanted. After all, the beginning and end of a rectangular pulse are switching processes that do not take place at one point in time but in a short time interval. It would be interesting to get information from Cornel about whether the linearity and inclination at the pulse edge are suitable models of the switching process. With a trapezoidal pulse, the inhomogeneity of the original equation is of course continuous and, in theory, the existence of the solution is then assured. However, since the expected inclination of the pulse edge is very close to the vertical, numerical problems are to be expected. The tangent will appear in the calculations. A sensitive divergence takes place near 90°. The length of the representable numbers will then be decisive.

The problem becomes even more interesting if equation solutions are generated using a sequence of trapezoidal pulses, with the pulse edges becoming increasingly steeper. Nothing can be said about the sequence of the associated equation solutions, whether it converges and thus trapezoidal impulses lead to acceptable approximate solutions. It would then be a very difficult problem to solve from real analysis (Lit.: Natanson, Hewit/Stromberg). The usual ideas of differential quotients and integrals would then no longer be sufficient.

Summary:
I still believe that with the tools available in MC only a "piece-by-piece solution" is possible. Regardless of this, I would be interested to know how the rectangular impulse is processed in the algorithm in the solution presented by Luc.

[Cornel]

I think in ttokoro picture he used V2(0s)=0V or something like a=-12 (having R1=10, C1=10^-8). Look at his picture and at my second picture from 2 comments below ttokoro's comment. But of course he should post better his worksheet.

I posted a second comment after his comment showing that with his chosen values for C1=10^-8 and R1=10, and having Adams/BDF for a=1 for example the V2(t) signal misses some pulses(more triangle in shape with this C1 value). I do not remember right now If for others values of parameter a, Adams/BDF misses pulses (I am on the phone now when writing this comment). And only for Radau V2(t) signal looks ok for all defined values of a=0...12 in this case of C1=10^-8 and R1=10. 

Anyway, for me is still unclear this behaviour of odesolve block, If is indeed as somehow Alfred comments or if Its a fault of odesolve block algorithm or has other reasons. That's why I suggested If someone can do this exercise at least in Mathcad 11 and Mathcad 15 for this moment to see If the results are the same as in Prime 10.