cancel
Showing results for 
Search instead for 
Did you mean: 
cancel
Showing results for 
Search instead for 
Did you mean: 

Community Tip - Need to share some code when posting a question or reply? Make sure to use the "Insert code sample" menu option. Learn more! X

Problems with solving a system of 5 differential equations 1st order

tsiebert
1-Newbie

Problems with solving a system of 5 differential equations 1st order

Hello PTC Community,

I have problems solving a system of 5 differential equations simultaneously by using ODESOLVE in Mathcad Prime 2.0. I tried other solvers too, but i have the most experience with odesolve. I just don't get the solver to work, no matter what I do. Could you please take a look at my file and tell me, if I have some Syntax-mistakes or anything like it.

The equations should be able to calculate the dispersion of gas, leaking out of a ruptured vessel in a cylindric form.

u represents the speed of the flowing gas

b is the radius of the cylinder

Theta is the angle between the direction of the gas and the horizon

c is the concentration of the leaking gas

T is the temperature of the flow

All those five dependent varibles depend on the independant variable sL, which represents the distance of the center of the cylindric flow to its origin at the ruptured vessel.

Am i maybe using the wrong solver? Any help or suggestions are greatly appreciated.

Tim

1 ACCEPTED SOLUTION

Accepted Solutions

In spite of my earlier comment I have had a go at rearranging the equations. The process I followed was:

1. Collect together all terms multiplying a given derivative on the RHS.

2. Define ancillary functions for each of these.

3. Because dtheta/dw doesn't depend on the other derivatives replace this in the other equations by the functions on its RHS.

4. There are three equations that now contain only du/dw, db/dw and dT/dw, so solve for each of these explicitly (ie solve for du/dw, db/dw and dT/dw)

5. Create functions from each of these solutions (they have a common denominator)

6. Now odesolve works.

See attached.

I will almost certainly have made a midtake somewhere along the line, so you will need to repeat the process, rather than relying on my results. I've plotted the results, but I've no idea if they look sensible or not.

Alan

PS I'll have a quick look at your latest attempts shortly.

View solution in original post

27 REPLIES 27

I couldn't spot a severe syntax error on first sight.

When I first open the file directly from you post I get an error telling me that the matrix of derivatives could not be created.

When I let the file recalculate or on saving it and reopening I get a different error: " Unknown error: exp%_too%_large".

So maybe its a bug in Prime or a limitation which should be dealt with. Prime is a work in progress, a rather unfinshed product and by no means can live up tp prior Mathcad versions like 15 or 14 (or 11 of course).

Unfortunately there is no way to save a file in MC15 format from within Prime and there is no converter available, so, as I'm not really that idle that I would retype the whole sheet, I cannot try if it would work in Mathcad 15.

Thank you for your fast answer,

yes those are the two errors i do get as well. Do you have any idea how i can find those limitations? Or do I have to look at each term in separate?

I am not sure. As you have possibly guessed I don't like Prime and so I lack experience using it. I'll stay with Mathcad 15 as long as possible.

I am not sure about the Prime bug or limitation and don't know how to verify unless someone from PTC will confirm.

Maybe you can try to simplify your huge expressions just to see if its a syntax issue. If the system gets somehow unsolveable at least the error message should change.

Its a pity we cannot try in MC15 as there is no converter available.

Hello Werner,

I was thinking, that maybe the solver starts at the right point, but it proceeds in the wrong direction. For Example if u take the speed u, it starts pretty fast at roundabout 400 m/s, but it will decrease now and not increase. Is there any way in mathcad, that i can tell the programme to go back and not forth, if you know what I mean.

Another way of solving this problem could be as u said the limitations. So if I want to do that, do I just define for example u(sL):=400m/s, 399.9m/s..3m/s before the solveblock, or is there a different way?

Regards,

Tim

AFAIK Mathcad takes the starting point from you initial conditions (in your case 0m) and the endpoint is set in the odesolve() command (in your file 100m). If the latter is smaller than the first it works backwards but it will always start at the points defined by your initial conditions. The solver cannot start at the "right point" as you demanded, as there is no "point" to start from - we don't know what u(100m) etc. will be.

Of course if you can define values for sL=100 you could change your initial conditions and use odesolve((),sL,0) to solve from 100 to 0.

Defining a range variable for sL should have no effect on the solving and would not make any sense.

And you cannot define a function as range as you sketched with u(sL):_....

My best guess (from the error messages thrown) still is that the problem has to do with the complexity of your equations. Not sure if MC15 would be able to deal with them.

So, I downloaded the trial version of Mathcad 15 and typed it in. While doing this, I found a big error in my fifth equation, so forget about my initial file.

I had to write the system without units, because odesolve is not able to deal with them. Furthermore, I had to addapt some indices, because as soon as i used the infinity-symbol, it produced an error, so my new index for surrounding is "1".

Problem is, I still got an error. It says: "This value must be a function, but has the form: any1."

Could you please have a look at this M15-file and tell me, if I made obvious mistakes?

Sorry, not a[1 but a.1 etc

You use in Prime a text index (a ctrl+_ 1), but in Mathcad 15 - a vector element a[1.

The text index in Mathcad 15 - a.1 (not a ctrl+_ 1 as in Prime).

Thank you for your advice,

I guess that was a pretty stupid mistake, but I blame my mathcad-Prime 2.0 experience (which is in some parts easier to handle by the way, a bit more intuitional) 😉

I corrected all indices, but without a result. Still the same error.

Divide by 0

0.png

Valery, what did you do to get Mathcad throwing that error?

My best guess is, that there is some kind of redundance, that the five equations are not independent.

The error message thrown was new to me in conjunction with odesolve. I saw the same message only in the form "....but has the form: unitless". So I tried to reproduce that very error with a much simpler system of equations and this is what I came up with

ODE_Error1.png

So maybe its a good idea to recheck the equations for consistency.

Its always the last vector element which is flagged red - it does not mean, that there is something wrong with that specific function.

You told us you had a mistake in one of your equations in your Prime sheet. When you corrected that, was the error message the same? The first message Prime throws, the one about not being able to create a table of derivatives, could correspond to an dependent system, I guess. The second, unknown error about an expression being too large lead me to think that it could have to do with the complexity of your system and an unknown intern Prime limit.

BTW, you can insert the infinity sign as index in a variable name in Mathcad 15, but its a little bit tricky:

1) insert the infinity sign anywhere on the sheet, select it and copy it to the clipboard (Ctrl-C)

2) Begin typing your definition, e.g. "u."

3) Press Ctrl-Shift-K and watch the cursor getting red. This is a special tricky mode which allows you to insert characters in a variable or function name, which are normally not allowed.

4) Now insert the infinity symbol from the clipboard (Ctrl-V)

5) Press Ctrl-Shift-K again to leave the special mode and keep on typing ":3"

uinf.png

Dispersion.PNG

Alan

Good Example, Alan.

In Prime2 my example gives an error about Prime needing more equations (so Prime's Odesolver is a bit cleverer than the one in MC15 according to more suitable error messages). Your example throws the same error which is seen the first time in Tim's file:

ode_error2.png

In contrary to Tim's solve block, recalculating the sheet does not change the error (In Tim's file we get an "unknown error - expression too large" then.

So it seems the problem with Tim's system is not due to an error in the ODE system but its simply that Mathcad and Prime are not capable enough to solve it. Is there a workaround? While I could create an example in MC15 where rearranging the terms made it work, I was not able to do so with your example. At least not with simple rearrangements. Substitution derivatives and rearrangement can do the job, but would be difficult with the 5 eq system. Next pic shows how sensible Mathcad is:

ode_error3.png

Werner Exinger wrote:

.... Is there a workaround?

The best approsch from Mathcad's point of view is to arrange the equations such that each equation contains only one first derivative. This can be done, in principle, for the dispersion problem because there are five linear equations in the five derivatives. Because each equation is large (they could do with tidying up by using seperate functions) this is a lot of work, and I don't intend to do it!

Alan

Hello Alan,

thank you, too, for your support so far. I also thought, that the derivatives on the right hand side would be a problem for mathcad, because in all the examples I found, there were none of that kind. What I then did, was that I tried to build up the system from one variable up to five. I made constants out of the variables and built up as many equations, as I could. When I couldn't add another equation, I tried to simplify that equation by leaving out terms from the numerator, until I could add it to the system.

I got to a system of 4 variables. Equation 1 is extremely simplified, equations 2-4 are complete. I left out equation 5 and I couldnt add it to the system, no matter what I did with it. Does this prove, that I can work with derivatives on the right hand side?

Regards,

Tim

P.S.: I found a typing mistake in equation 1, I added a new corrected file (It did not change the result of receiving an error).

Does this prove, that I can work with derivatives on the right hand side?

It is possible to have derivatives on the RHS (see Alans first example), but Mathcad seems not to like them, as Alan wrote, and ever so often will refuse to do its job. The exact conditions when MC will accepts that kind of equations and when it fails are unclear to me.

Find attached a file which follows Alans idea of solving your equations for the derivatives first. Of course I haven't done it by hand but let MC do the work.

I included Alans example in front of the file, which failed with the same error as yours and the method worked fine for it.

When I try to apply it to your system, Mathcad crashed on my machine. It could be because of lack of memory, I could imagine. Also the solving for the derivatives is very slow, I haven't timed it but I guess it took about an hour on my machine (which is rather slow and has only 1 G of mem, so you may have more luck).

I also ran into a strange error of Mathcad not being willing to make substitutions with a function name consisting of a Greel character - fortunately there is a workaround (other than renaming the function):

GreekSubst1.png

Wow, thank you so much Werner. I really appreciate your work.

I will immediately have a look at your file. Hope this works, I will keep you updated;-)

In spite of my earlier comment I have had a go at rearranging the equations. The process I followed was:

1. Collect together all terms multiplying a given derivative on the RHS.

2. Define ancillary functions for each of these.

3. Because dtheta/dw doesn't depend on the other derivatives replace this in the other equations by the functions on its RHS.

4. There are three equations that now contain only du/dw, db/dw and dT/dw, so solve for each of these explicitly (ie solve for du/dw, db/dw and dT/dw)

5. Create functions from each of these solutions (they have a common denominator)

6. Now odesolve works.

See attached.

I will almost certainly have made a midtake somewhere along the line, so you will need to repeat the process, rather than relying on my results. I've plotted the results, but I've no idea if they look sensible or not.

Alan

PS I'll have a quick look at your latest attempts shortly.

Hello Alan,

So I rewrote most of your file myself and i have to say I really appreciate your help here, you've put quite some work into this and I wanted to thank you again. I still have 2 questions, though:

- How can you leave out the independant variable w in most of the terms and in the end you can just put them into the variables again? How is this possible? I see the problem, that when you define something in dependancy of w, there is an error saying "not defined Variable", but I can't explain why it works out fine in the end.

- How do you handle dT/dw, because in the formulas, there is no dT/dw, but d(Roh)/dw. I get stuck on part 5 of your way through solving the equations. Do you just multiply T/T and end up with: dT/dw*(Roh/T) ? I thought that you can't do that, but that would be my only explanation. I didn't try everything in my mind yet, but this is the point, where I cant comprehend your way.

Tim

Edit: All right, I understand, what you did with d(Roh)/dw. It's the same as p/R*d(1/T)/dw. This is then equal to -p1/(R1*T^2)*dT/dw! I think then you re-substituted it to -Roh/T*dT/dw.

Tim Siebert wrote:

.... Does this prove, that I can work with derivatives on the right hand side?

It certainly seems so, though because you've set T(w) to a constant, the derivative terms involving T(w) will go to zero. However, Mathcad seems pretty picky about when it will and when it won't accept this situation, so I still recommend having no more than one equation with derivatives on the RHS.

Incidentally, in your equation for db/dw there are two terms involving du/dw. In one it is multiplied by a10.a11.p1/R1T(w). In the other it is multiplied by a10+a11.p1/R1T(w). In my worksheet I assumed these should both look the same. If I'm wrong about this that's one error that must be corrected in my worksheet!

Alan

Hey Alan,

thank you again for your support, you are a really big help! So I looked it up and you are right, that is a typing-mistake. It should both be a10+a11.p1/R1T(w). I will now need some time to comprehend your work. But I do think that this is the right path.

Regards,

Tim

But I do think that this is the right path.

I, too, think that Alans sheet could be the solution to your problem.

I just tried mine again, this time all went OK (solving took just 8 minutes) but Odesolve throws an error and I have no clue whats wrong this time.

Interesting point, as Werner, I really would like to know how you produced that error.

What can I do to prevent this? Is there a way to limit parts of the functions?

We have problems with odesolve:

- in Mathcad 15 we can not work with units

- in Prime we have problems with lables of variables, functions and constants (find and see please my messages here about it).

And some litle remarks

- g is built-in constant

- do not use the = operator after the odesolve function

- in Mathcad 15 we can not work with units

which is rather annoying indeed as we have to use some nasty and unsatisfyable workarounds.

- in Prime we have problems with lables of variables, functions and constants (find and see please my messages here about it).

I remember that problem from some other threads. I looked over the expressions (but not with the necessary meticulousness, though) and didn't spot anything of that kind. Did you find some of these inconsistancies in Tim's file?

According to the error mesages thrown I guess the error to be something else.

- g is built-in constant

Tim redefined it at the top which should not do any harm.

- do not use the = operator after the odesolve function

Agreed, but deleting it does not change the situation.

Do you have the same problem in your Prime3 installation?

Thank you, too for your fast answer.

Would you recommend to use a different type of solver, if odesolve has some known "bugs"?

I will look for your previous posting about the problems, thank you so far!

Sometime I use an open program to see a part of ODEs solution:

http://twt.mpei.ac.ru/ochkov/Mathcad_12/5_13_My_rkfixed.png

Hello Valery,

I am sorry, but I do not understand, what you are trying to tell me. Is this just an example of your open program, or is it specifically a part of one of my equation?

So your goal here is to find parts of an equation, that tend to go to infinity, am I right?

Top Tags