I played around a bit with the examples by Radovan.

My observations:
- using stack() or not is irrelevant in the single first order equation case.
- D and initial condition may be scalar or vector-valued.
- if() and Max() and eval() work, if D is defined inline (in the solver argument list), they do not work, if D is defined outside.
- the different behaviour of Rkadapt is due to poor convergence of the adaptivity algorithm.
- if rkfixed (and possibly the other solvers as well) fails, then you can evaluate the expression using = and get an error
- if rkfixed succeeds, there is no way to evaluate it using = (key is ignored), you have to use assignment instead

---

if-max-eval.sm (12 KiB) downloaded 117 time(s).

Updated. Source code refactoring.

The problem here is not the solvers, the problem is in the values ​​of the function parameters.

[albumimg]263[/albumimg]

Wrote

I played around a bit with the examples by Radovan.

My observations:
- using stack() or not is irrelevant in the single first order equation case.
- D and initial condition may be scalar or vector-valued.
...

Do not recommend using scalars. The point is that the ode functions in SMath Studio are described as follows:

Quote

new ArgumentInfo[] {
new ArgumentInfo( ArgumentSections.ColumnVector ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.Function )
}
);

For now there is no strict monitoring parameter types, but it can occur in the future. I can make support for scalar parameters, but then for every function I need to write it analog with different set of types of parameters.

Wrote

Do not recommend using scalars. The point is that the ode functions in SMath Studio are described as follows:

Quote

new ArgumentInfo[] {
new ArgumentInfo( ArgumentSections.ColumnVector ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.RealNumber ),
new ArgumentInfo( ArgumentSections.Function )
}
);

For now there is no strict monitoring parameter types, but it can occur in the future. I can make support for scalar parameters, but then for every function I need to write it analog with different set of types of parameters.

Ok, I see. How about a generic ODESolve(method, y0, t0, t1, n, D), where method is some optional string indicating the numerical procedure. Then, whatever convenience UI stuff you add goes to all solvers.

Similar to the FindRoot function you might adopt the embedded assignment option, e.g. by specifying the initial conditions and iniial value of the independent variable in equation form like

[MATH lang=eng]ODESolve(mat(v≡1,s≡0,2,3,1),t≡0,1,10,D(t,y))[/MATH]

That would create vectors v, s and t, ready for plotting or further processing. In case of two requested vectors you might augment them into something like v_s which would immediately be good for plotting.

Dear uni and omorr

I can not express in words how thankful I am and my admiration for your work. SMath studio will be my preferred tool from now on.

To help others I have the intention of making a “stirred tank reactor” example with microorganisms growing and consuming a substrate; I already have such an example in Mathcad from the time I was teaching. There is already a similar example posted in https://smath.com/wiki/Examples.ashx, but an example with the new plugin functions would be useful I believe.

Thanks again!

You are welcome and thank you
I am looking forward to your examples.

Regards,
Radovan

I tried to summarize the results of the recent discussions as to how to use the ODE solvers. Up to now I assume that everything might work with any of the solvers but that is not tested thoroughly. It would be desirable to have some sort of benchmark problems and to compare the performance of the individual solvers. That might result in sort of usage recommendations, if not in a generic solver just like FindRoot.

Feel free to challenge and correct the statements in the attached handbook pages.

Edit: Minor corrections in the introduction (example description had a sign error) and in the requirements for D (there was a spurious result display for an expression with embedded :=, which I could not reproduce)

---

section math ode.zip (66 KiB) downloaded 131 time(s).

Martin,

Once more, all my respect to you for having a good will and patience to extract all these things and put them together. I have to admit that sometimes I envy you for that

Regards,
Radovan

Updated. Added functions from the Intel® Ordinary Differential Equations Solver Library: rkm9st(), mk52lfn(), mk52lfa(), rkm9mkn(), rkm9mka().

mk52lfa() not finished yet.

[albumimg]272[/albumimg]

---

ODESolvers. Van der Pol oscillator.sm (10 KiB) downloaded 148 time(s).

I have never imagined that so many ODE solvers would be available in SMath. I wish you to stay inspired because you are keeping our positive spirit high
Many thanks for that.

By the way, are you going to to introduce some additional parameters like RelTol,AbsTol in DotNumerics. As I could see from the Reference Manuel for these solvers, among many other parameters (I suppose you set them by default values) there are parameters like ep, tr for absolute and relative tolerance.

Regards,
Radovan

Wrote

I have never imagined that so many ODE solvers would be available in SMath. I wish you to stay inspired because you are keeping our positive spirit high
Many thanks for that.

By the way, are you going to to introduce some additional parameters like RelTol,AbsTol in DotNumerics. As I could see from the Reference Manuel for these solvers, among many other parameters (I suppose you set them by default values) there are parameters like ep, tr for absolute and relative tolerance.

Regards,
Radovan

You can use RelTol and AbsTol in the same way. I included them in to the plugin. This code is running in the unmanaged space. Therefore, it is platform-specific. During loading the plugin from its resources on disk is copied library for specific platform (iode.dll). I'm wondering to know does it work in 64-bit system?

Now I can use any function that is written in C/C++ or other languages.

Works on x64 systems (at least on my notebook), except for the mk52lfa(Y,t0,tmax,n,D(t,y,4)) (Y not defined)


best regards,

Davide

Wrote

Works on x64 systems (at least on my notebook), except for the mk52lfa(Y,t0,tmax,n,D(t,y,4)) (Y not defined)

best regards,

Davide

Thank you, Davide. Unfortunately, this function requires the Jacobian. It must be specified explicitly. Some other functions also require it, but it can be omitted. I'll have to change the number of parameters to a function call.

Wrote

You can use RelTol and AbsTol in the same way. I included them in to the plugin. This code is running in the unmanaged space. Therefore, it is platform-specific. During loading the plugin from its resources on disk is copied library for specific platform (iode.dll). I'm wondering to know does it work in 64-bit system?

You are right. RelTol and AbsTol worked with these functions as well . I tried it and can confirm that mk52lfn() and rkm9mkn() are quite fast and worked as expected with Abstol=Reltol=10^-7 . There is no secret which one of the examples i tried . This worked on my 64-bit laptop. I also could see the iode.dll in the plugins folder. I used few SMath SMath versions (including few ones by Martin) and am not sure if 32-bit is enabled or not. Anyway, I think it is not so important at the moment

Wrote

Now I can use any function that is written in C/C++ or other languages.

All my respect for that . I hope that you will find a nice "candy store" now and allow us to taste it as well

Regards,
Radovan

Here is a nice little toy for exploring precision control. Have fun!

---

Wrote

Unfortunately, this function requires the Jacobian. It must be specified explicitly. Some other functions also require it, but it can be omitted. I'll have to change the number of parameters to a function call.

Does that mean that currently, the routines use an internal approximation of the Jacobi matrix, whenever they need it? Could the matrix be generated analytically or numerically by the functions from Davide's Nonlinear Solvers? Would that require user input?

Dear colleagues
I made an example on how SMath ODE solvers may be used.
Comments and improvements are welcome
If you thing it would be useful, feel free to post it or use it wherever it is best suited.
I am not sure of what would be appropriate.

---

Box_models.sm (108 KiB) downloaded 156 time(s).

Nice example

As a chemical engineer I am quite "attached" to chemical, biochemical reactors . That is one of the reason I was nagging before about ode solvers (at the moment we have many of them - thanking to uni .
Quite a chance that this model will fail with some nonstiff solvers. Fortunately, and thanking to uni, there are three plugins (at the moment) with ode solvers (OdeSolvers, DotNumerics and Mathlab C++ Math library). You can try all of them.

Regards,
Radovan

I have solved one equation in sage math but not able understand how to solve in smath anybody can guide me on this. below is my equation


exercise1

r = var('r')
PI = pi.n()
rc = 281
rw = 262
Bc = 36/180*PI
F = (r^2-rw^2)*sqrt(rc^2-r^2*sin(Bc)^2)/sqrt((rc^2-rw^2)^2*rc^2*cos(Bc)^2-(r^2-rw^2)^2*(rc^2-r^2*sin(Bc)^2))
F.show()




z = function('z',r)
d = desolve_rk4(diff(z,r)-F,z,ics=[rc,0],end_points=255,step=-0.01)
list_plot(d,plotjoined=True, aspect_ratio=1)

---

From the start, you are all right and all wrong !

1. You have to understand what you are doing and how it can be done numerically.
2. No need for all these Moulton gadget, Fhelberg ... etc. For instance: Fhelberg
extends the RK domain where it has no meaning [Trust my years old testing].
3. All those librairies, proven ? only few of them [NAG, ACM... OK ? !]
4. My Smath version came with Maxima plugins: some work some don't. On that, Smath
is not an integrated version. Your clients will be soon fed up of shoping for zillions
of unproven plugins and eventual incompatibility between versions.
5. On the other hand, as long as Smath is 32 bits it is not really a CAS.
6. Too much automation in solvers may exclude finding a solution. What can be done
with little manual pushing may take days to de-automate. In Mathcad 11, Levenberg-Marquardt
works generally well but many fits have to start manually, then Newton, then CG, refined LM.
6. All those 3D plots from "Uni" make me envious, but why should I scratch my right ear with
my left hand or ask someone else to do it [unless paralysed]. The GNUplot works, but that's
the only one. Martin couldn't put color in my famous "Breather".

Very interesting.

Cheers, Jean