Wrote

Hmm... why dirty? Do not understand, sorry. This looks quite regular to me, except that it is not flexible. The third parameter must always have the same name as in the definition of D().

Regards,
Radovan

The method is not reliable. If the parameters have values in advance (before the definition of D(x,t)), then the approach works if the definition is set to numeric optimization and does not work if it is set to symbolic optimization. As you can see in the images, symbolic definition leads to replacement of parameters by their values. This does not happen with numeric optimization.

Would you have guessed that? I don't like methods, which make me doubt whether I know what I do in SMath.
[albumimg]146[/albumimg][albumimg]147[/albumimg]

I see. To be honest, I might guess you were thinking on this quite problematic, and sometimes frustrating behavior.

This was already mentioned more than few time, here and there on the Forum. The complaints ranged from different results to rather different calculation times etc. This numerical-symbolic conflict in SMath must be resolved in one way or another.

Regards,
Radovan

Hello uni,

my impression is that Rkadapt and rkfixed both give equidistant time steps in the output. Is there a way to access the internal adapted timesteps (just for demonstration)? I'd like to plot the time step vs time in order to show the adaptivity, but failed to do so up to now.

Martin

[albumimg]149[/albumimg]

---

Lorentz-Atttraktor.sm (16 KiB) downloaded 136 time(s).

No, there is no way. The original algorithm has an adaptive step, but I wanted to do it by analogy with similar functions in Mathcad (rkfixed and Rkadapt), where they both give the result on a fixed grid. It took a while, but I figured out how to make an adaptive algorithm to do the calculation in the fixed points. To test the results, there is only one possibility - to compare the accuracy of integration (rkfixed and Rkadapt). For my tasks was no need for adaptive step.

Could you take a look at this library, LibSBMLSim ?
It could be interesting for the ODESolver plugin. It is LGPL licensed and provides a C# binding.
http://fun.bio.keio.ac.jp/software/libsbmlsim/

Edit:

Two more libraries with odepack implemented:

SmartMathLibrary.OdePack.dll

SciPy for .NET via IronPython
I think this approach needs .NET4.0

Edit:

NumPy and SciPy for IronPython / .Net

Dear Viacheslav
I tried your plugin, and I have compared the results with the equivalent function in MathCAD, see attached file.
Either something is not working or I have misunderstood how the function should be called.

I downloaded the ODEsolvers from
ftp://82.193.156.30/SMath/0.96.4909/plugins/

with the portable SMath Studio 0.96 (build: 4909).

Could you please tell me what is wrong?
Thanks in advance.
Best regards Ignasi

---

Test_Rkadapt.sm (16 KiB) downloaded 123 time(s).

Hello Ignasi,

These kinetic equations of this kind are often troublesome and this is also a "tricky" one. In these cases a stiff solver might get the problem resolved - Rkadapt is not a stiff solver.
The problem causes the component B because its concentration tends to be small (k2,k3 >> k1). Therefore, it is advisable to try to solve these kind of problems first without the intervention about negative concentrations. I tried this and the file is attached. On the other hand, the same problem was solved with Mathcad. As expected, Rkadapt in Mathcad performed much better. Rkadapt in Smath needed rather high integration steps for a bit larger reaction time, but it could not resolve the problem with larger reaction time as Rkadapt-Mathcad did. Rkadapt in Mathcad failed here too for a larger reaction time, and I used the stiff solvers - see attached.
[albumimg]227[/albumimg][albumimg]226[/albumimg]

Regards,
Radovan

---

Test_Rkadapt-corr.sm (12 KiB) downloaded 124 time(s).

Thank you, Radovan. I think you're right about the different concentrations. The algorithm tries to find step to achieve accuracy in the integration, but it is impossible because of the different absolute values ​​of derivatives. Unfortunately, I don't quite understand the chemical kinetics and can't explain the results. The examples that I've seen before on this subject have been much easier.

EDIT: It seems I found sources for stiff ode solvers: DotNumerics.

Thank you uni for the reply and for the consideration to take into account some stiff ode solver

Those stiff situation are quite often in chemical kinetics. It seems that this is actually the common situation. By the way, as I mentioned few times before, any stiff ode solver is quite appreciated in SMath. If you decide to apply the mentioned ones, you can use this example as a test one. Moreover, I often tested many solvers with this example Bioreactions. Ode solver should resolve this and give the concentration profile as in the next post but without the intervention on zero concentration.

Regards,
Radovan

Wrote

Dear Viacheslav
I tried your plugin, and I have compared the results with the equivalent function in MathCAD, see attached file.
Either something is not working or I have misunderstood how the function should be called.

I downloaded the ODEsolvers from
ftp://82.193.156.30/SMath/0.96.4909/plugins/

with the portable SMath Studio 0.96 (build: 4909).

Could you please tell me what is wrong?
Thanks in advance.
Best regards Ignasi

Try now with the DotNumerics plugin from the Gallery.

EDIT: Ignasi, your example is used in the DotNumerics to illustrate how the function works.

Quote

// Problem from chemical kinetics. Consists of the following three rate equations:
// dy1/dt = -.04*y1 + 1.e4*y2*y3
// dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
// dy3/dt = 3.e7*y2**2
// on the interval from t = 0.0 to t = 4.e10, with initial conditions
// y1 = 1.0, y2 = y3 = 0. The problem is stiff.

[albumimg]232[/albumimg] [albumimg]239[/albumimg] [albumimg]240[/albumimg] [albumimg]241[/albumimg] [albumimg]242[/albumimg]

---

GearsBDF. Example 2.sm (9 KiB) downloaded 125 time(s).

Thank you uni very much. I was pleasantly surprised how fast the dn_GearsBDF() was comparing to the other two solvers Rkadapt() and al_rckrkadapt(), and the success regarding the problematic kinetic equations given by ignasi. I can also see that my mentioned problem was not solved because it goes to negative y2 results. It is a problematic point because of the close to 0/0 situation. It is not that problematic. As I said, I could see some other software included stiff solvers which can deal with that (implemented some kind of Gear algorithm or lsode based), but most of them gave the same result as presented. I really do not know what was the catch about this. Moreover, I just tested the three solvers on Lorentz attractor equations. It seems that dn_GearsBDF() gave some different result for the second picture - see attached. Do not know the reason.
[albumimg]238[/albumimg]

Regards,
Radovan

---

ODE-Lorentz-dn.sm (11 KiB) downloaded 105 time(s).

I think this is normal. It caught second petal. Perhaps the first two cases are not correct. As for Ignasi example, it is taken from the book: Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems.

dn_GearsBDF() not finished yet. There are some little issues.

---

Thank you uni,

dn_GearsBDF might be the ODE solver of choice for SMath. I would not mind some more solvers as well

Regards,
Radovan

Dear uni:
First of all, you deserve many thanks and admiration for making available these plugins! They are exactly what many scientists need. I tried your new functions as soon as I read your post, naturally.

By the way, you are right about “my example”: it is a classic example found in several places. I found it in the following document: http://www.ntis.gov/search/product.aspx?ABBR=N8024405
Surprisingly, Rkadapt in my Mathcad could cope with it.

Nevertheless I still got a problem. It is related to rates of processes that only may go in one direction. For example, if you have gasoline and oxygen from air your car will move forward and produce CO2 and H2O, but if you start with CO2 and H2O the process will not go in the other direction. Many biogeochemical processes catalysed by microbes are of such kind. In order to simulate these systems one could use the “Max” function, or an “if” statement.

But the examples in the attached file do not work.

Two years ago omorr wrote (#9)

Quote

I tried to define function D as a multiline function but this will not work”, and I believe that the following posts up to #14 were about this subject. As far as I understand the comment from omorr in #14 is still valid? “It seems, as you mentioned, that Rkadapt will not work with multiline function. . . . I also tried to define function f(x,t) by using "eval" - it does not work either with rkfixed or Rkadapt. Due to my previous experience, I thought it would not work. I know that this kind of functions, like you did with rkfixed and Rkadapt, will have troubles due to the SMath way of working

I feel somewhat ashamed to ask the following questions taking up some of your precious time, but I would be very thankful if you could help !

Questions:

1. Is it expected that “Max” would not work? Or is it a problem easily fixed?
   
2. Can the function be defined as a multi-line or not?
   
3. In your post #15 you write “Updated description.” I tried to figure out what did you mean, and search for a description. Where is it?
   
4. Is it still true that:

   Quote

   ...the side effect is that rkfixed and Rkadapt will be working wrong if this function (D(t,y)) uses t and vector y with t,y previously defined” (omorr, #9)?

Thanks in advance! Ignasi

---

Max_and_if.sm (10 KiB) downloaded 109 time(s).

Hello Ignasi,

I hope you would not mind if I try to give you some answers (I hope I am right about all of them).
By the way, your example is with two equations, but you actually need one (you made some typos as well - mixed lowercase y, upercase Y variables). The derivative of the second is zero. This was a trick in Mathcad if you wanted to use some parameters. In Smath you do not need that - see the examples in this thread please.

Wrote

Questions:

1. Is it expected that “Max” would not work? Or is it a problem easily fixed?
   

See the picture please. As far as I know, the comparison thing will work as explained, otherwise will not. I really do not know if there could be a fix inside SMath in order to use Max(), if, or eval(). There is a continuing problem regarding this numeric-symbolic issues and I do not know if there would these things ever be solved.
[albumimg]260[/albumimg]

Wrote

Can the function be defined as a multi-line or not?

Yes, it can. See the examples and pictures please

Wrote

Is it still true that: ...the side effect is that rkfixed and Rkadapt will be working wrong if this function (D(t,y)) uses t and vector y with t,y previously defined” (omorr, #9)?

Yes, it is. The first two parameters of the D() function should be symbolic and not defined previously - check them first.

Regards,
Radovan

---

Max_and_if_corr.sm (8 KiB) downloaded 109 time(s).

Just to mention one more way that works: the internal function max() is int()- and ode solver-proof:
[MATH lang=eng]D(t,y):max(mat(nthroot(el(y,1),a),0,1,2))[/MATH]

Wrote

Just to mention one more way that works: the internal function max() is int()- and ode solver-proof:
[MATH lang=eng]D(t,y):max(mat(nthroot(el(y,1),a),0,1,2))[/MATH]

Tnank you Martin, . I forgot about it. As Davide Carpi once said - everything could be done

Thank you omorr and mkraska
I noticed too late of my misprints (I changed my file at the lust minute without double checking). And yes, I understand that I did not need two variables. In the system I am working with I have six variables, and I tried to scale down my system to small and clear examples to illustrate the problem I was encountering. I will double check next time.

Thank you very much omorr for pointing out a solution by using a logical condition!! That is elegant. Also thank you very much for your explanations and confirming how the plugins are expected to work.

The question of the multi-line with the plugins: rkfixed, Rkadapt, and dn_*.

1- I have looked at the examples in this topic, and first of all I have to confess my ignorance, for example, omorr, in your post #14 you define f(x,t) with a program line. Also D(t,z) in post #30 in file "Test_GearsBDF". I do not understand the purpose of these lines, and I have not been able to find an explanation in tutorials or help files. Is this explained somewhere? In the case of "Test_GearsBDF" the line may be removed without any effect. Can you confirm that the purpose of the line was just to show and test that a multi-line definition was possible at all?

2- I have tried your example in #14, and it does not work with my portable SMath 0.96.4909 (as expected I assume) neither with rkfixed or with Rkadapt. The file "Test_GearsBDF" posted by uni in #30 works. My conclusion is then that a multi-line is possible, but an "if_else_" is not possible. Please let me know if I am wrong.

Radovan said almost everything. Add on my own, the solver works with matrices only for now. You must use the stack() if there is no system of equations. The same is true for the initial conditions. As for examples, they work fine with other conditions.

EDIT: Added example with multiline D(t,y). This is a real example which illustrates the practical use of the D() function as a program unit.

[albumimg]261[/albumimg]

---

DorNumerics. ODE. Arenstorf orbit.sm (91 KiB) downloaded 125 time(s).

Wrote

The question of the multi-line with the plugins: rkfixed, Rkadapt, and dn_*.

1- I have looked at the examples in this topic, and first of all I have to confess my ignorance, for example, omorr, in your post #14 you define f(x,t) with a program line. Also D(t,z) in post #30 in file "Test_GearsBDF". I do not understand the purpose of these lines, and I have not been able to find an explanation in tutorials or help files. Is this explained somewhere? In the case of "Test_GearsBDF" the line may be removed without any effect. Can you confirm that the purpose of the line was just to show and test that a multi-line definition was possible at all?

2- I have tried your example in #14, and it does not work with my portable SMath 0.96.4909 (as expected I assume) neither with rkfixed or with Rkadapt. The file "Test_GearsBDF" posted by uni in #30 works. My conclusion is then that a multi-line is possible, but an "if_else_" is not possible. Please let me know if I am wrong.

I just tried these things to test it and to figure out how to solve some problems. This is the case about multiline functions that you asked. They are not needed sometimes. You can do the same thing with or without them. To be honest, I still am not sure what is going to work and what is not. I am just trying and sometimes it is working, and sometimes not - quite often a trial and error procedure. I am sorry, but I can not say why is this happening. For instance, uni explained about that the initial conditions and D() function must be in vector form (returning vector) even it has only one returning value. He mentioned stack() function to accomplish this. I know that stack() with even a single scalar element will return a 1x1 vector, and to be honest, I was surprised (again) that these things with Max(),if and even eval() will work by using stack(). See the picture please, by using your example.
[albumimg]262[/albumimg]
I can not explain why stack() inside dn_*() and rkfixed() will work but outside will not work. I can not explain why Rkadapt() fails here. Anyway, I have no other choice but to try and to see what is going to be. I hope that all these things will be find and explained in one place (maybe in the Help file made by Martin ). At the moment, all the user experiences are scattered on this Forum.

Regards,
Radovan

---

Max_and_if_corr_2.sm (17 KiB) downloaded 119 time(s).