My first mathcad worksheet- many glitches

[KP_10752086]

I am using Mathcad Prime 5.0.0.0, Windows 11, for the first time as instructor in a civil engineering course having taken over from the previous instructor mid-semester.  (Ordinarily I use Matlab.)  My worksheet is supposed to solve the differential equation for liquid spilling out of a small hole in a tank modeled by the ODE dh/dt=-K*sqrt(h) in two ways: odesolve and rkfixed.  Here are the glitches I encountered:

1. Plotting functions h(t) (the exact closed-form solution) and q(t) (the output from odesolve) did not work for me.  I received an error message about complex or infinite values.  However, it is clear from its formula that h(t) cannot take on infinite or complex values.  In addition, Mathcad repeatedly changed the units of the t axis from s (which I typed) to 1/s * s, which is incorrect.

2. The exact solution to the problem exists and is well defined up to about t=85000,  After that point, the tank is empty and the ODE is undefined.  However, odesolve was able to integrate only to about t=60000.  See the glitch with z(t) in the attached worksheet.

3,  I also solved the IVP with rkfixed.  The units came out incorrectly.  The dependent variable column of the output was labeled m*s instead of m.

I'd appreciate any help since I need to hand out this assignment soon.  I've attached my worksheet.

[LucMeekes]

1. You did not define t. Hence Prime uses its automatic definition, which is from -10 to 10, without units.

If you define t as (e.g.):

before the plot, then it plots as expected.

Success!

Luc

[KP_10752086]

Great, thanks! This fixed obtaining plots. Do you know the answer to my other questions as well?

[LucMeekes]

No.

In (real) Mathcad 11 the odesolve runs just fine:

I've checked Prime 7 and 8, they both quit at (about):

where

I have no clue as to why.

Maybe someone with Prime 9 can check....

Success!
Luc

[KP_10752086]

Thank you! It looks like something good was done in v. 11. How about #3? Why do I get the wrong units? Thanks again for the fast responses.

[Werner_E · Accepted solution]

---

@KP_10752086 wrote:

Thank you! It looks like something good was done in v. 11.

That for sure. The old Mathcad 11 is considered by many users as being the best Mathcad version ever built. This opinion is based to a good part on the fact, that in this version there still is integrated a slimmed down version of Maple for the symbolic operations.

I have tried your odesolve block in Prime without units with no success and the same ODE also fails in Mathcad 15 with the error message "This function cannot be used here" without indicating which function this text refers to.
BTW, of course the solve block fails in the current Prime 9, too - so the problem still is not fixed.

[Werner_E]

Your first question is already answered by Luc. You can't use the quickplot feature (plotting without defining the abscissa variable) with units - you have to define a range with units for the abscissa values.

I would add, that you may consider formatting the abscissa as hours instead of seconds. You should also know that on each axis you can edit the first, second and last scale value (can be quite a fiddling around). This is often necessary as Primes auto-scaling not exactly the best you can imagine.

I am not sure but the OdeSolve() block seems to fail because of numerical round-offs / inaccuracies. It looks like the numerical value of z (or q) is getting negative and so its the root which actually throws this error.
You may consider one of the two workarounds you can see in the picture:

The differences can be seen in the next pic and I guess its in a region which is not of interest to you anyway

In the workaround the "max" function could also be replaced by q*(t)=-K*Re(sqrt(q(t)) with the same effect.

The wrong unit returned by rkfixed is strange.

According to the help (the example there) the initial value should be a vector (a 1x1 vector in your case).
Providing this 1x1  vector instead of a simple scalar cures the wrong units problem (have no idea as to why it does, though).

The same example in the help
-> Example: Solving a First-Order ODE Initial Value Problem

also treats the second argument of the vector function as a vector even in case of a single function first order ODE

But this seems not to matter as your sheet is working both ways

   or    

Worksheet in format Prime 5 attached

[KP_10752086]

Interesting. I knew the solution might become negative due to roundoff near the end, which is why I ended rkfixed at tend-1000s. I guess I don't know what odesolve does. Is it symbolic?

Thank you.

[KP_10752086]

Ηi. Can you please tell me what odesolve does exactly? Should I even be teaching rkfixed? I wish the documentation offered more information on the methods.

Many thanks. 

[Werner_E]

Odesolve is pure numeric.Its similar to AdamsBDF, but it does not return vectors but a function (which internally is just interpolation the vector values).

And yes, more information about the numeric algorithms under the hood would be welcome. In real Mathcad we had the option to chose among a variety of algorithms odesolve() would use - default was Adams/BDF, others to chose were Fixed, Adaptive and Radau. Some numeric algorithms used for solve block were changed, when PTC began burying the real Mathcad and putting the weaker Prime in its place and as some examples here in the forum already showed, that wasn't always an improvement. There is no comprehensive documentation of the used algorithms as far as I am aware. The Prime help only states that OdeSolve uses Adams/BDF algorithm - the hybrid algorithm which can be used for non-stiff as well as stiff ODEs.

BTW, Mathcad unfortunately never was able to solve ODEs symbolically out of the box. Starting with Prime 9 PTC has added the functionality to symbolically solve some basic ODEs and promised to improve that feature.
See here:-> Solving ODEs Symbolically

Basically your ODE should be simple enough for P9 to solve symbolically, but it fails because of the usage of units - the symbolics does not know anything about units

After getting rid of the units the result is not really useful, either

and here is what we get using no numbers at all

Seems like PTC has a lot of work to do and this isn't the only thing that needs fixing...