solve() and roots() - again :(

solve() and roots() - again :( - Messages

#1 Posted: 4/8/2012 5:06:30 PM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Hello,

Unfortunately, there are still continuing problems with the nonlinear solvers.
Look at this example please Primer65a

It is a problem with moderate complexity. The goal is to solve a nonlinear equation.
Here is a screenshot from the file - v0.93


As it could be seen solve() and roots() gave similar results (T~328) (green and blue colored regions) but look at the function values - quite far from zero(~7*10^5). Something must be wrong with the convergence criteria of these functions, or I made some mistake I could not see.

It was enough here just to make few secant iterations and to get the root (red colored regions) (T~343) with the function value quite closer to zero (~10^-7).

I must admit that it took me quite a lot of time to get this result. Calculations took quite a long time (sometimes more than 20min) - different optimization did not help - until I realized that eval() on the right place will dramatically speed up the calculation (just few seconds).

I really do not know why eval() in that procedure made this calculation much, much faster. It seems that eval() will help in speeding the calculation but I did not realized yet when to use it. It is a bit "trial and error" procedure.

I commented this few times before, and I hope you would not mind if I mention this again.
Efficient nonlinear solvers (and some other numerical methods) are essential and crucial methods for every math software. It is a "must have" feature. Thanking to the Andrey's efforts there are solve() and roots() in SMath. Andrey did his best to make them but, unfortunately, in spite of his efforts it seems they are not god enough. I am not a programmer, have a vague idea how SMath plugins are working, have no idea how Andrey made solve() and roots() - I would not understand it anyway. I just feel (following some other math software) that the only solution is incorporating some established numerical libraries (made in C,C++,C# etc.) into SMath like "black boxes". For instance, I used Fortran and C procedures without knowing how are they working - just read what are the input parameters and what they return as results. I am just hoping that SMath plugin system can provide something like this and someone will do that eventually.

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#2 Posted: 4/10/2012 2:39:08 PM
kilele

kilele

133 likes in 397 posts.

Group: User

Radovan, If I might venture an opinion, you seem to be the appropriate person to help Andrey by writing some algorithms in the form of PSEUDOCODE for solving nonlinear systems
#3 Posted: 4/10/2012 3:31:27 PM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Hello,
Wrote

Radovan, If I might venture an opinion, you seem to be the appropriate person to help Andrey by writing some algorithms in the form of PSEUDOCODE for solving nonlinear systems


I hope you would not mind, but I have to disagree. There is already lots of libraries and code for numerical calculation out there and freely available. There are mature, well established, optimized and thoroughly tested numerical procedures made in different programming languages. It would be kind of wasting the time if we do it again, like "inventing a wheel". Of course, we can use that code in its original form but imagine using some well known Fortran or C++ code for, say, solving different kind of nonlinear or differential equations. One of the logical reasons for using those compiled languages is when you need power and speed, when you have huge amount of data, thousands of equations and quite complex problems. For some problems of moderate complexity using Fortran of C++ is hardly justified. You will spend much of your time on programming issues than on your problem solving. Therefore, in those cases we just need some attractive math environment like SMath equipped with lots of numerical procedures (made from the above mentioned libraries via the plugin mechanism). I do not understand how the plugin mechanism is working but looking at the environments which can be used for making them IDEs and programming languages for developing plugins I suppose that the SMath plugins for numerical calculations can be made by using some of those languages and already made numerical libraries in the same languages.

Of course, as I am not an expert in these things I might be totally wrong and this might not go that easy. Some plugin makers can give us more accurate and precise answers on this issue.

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#4 Posted: 4/10/2012 6:47:30 PM
kilele

kilele

133 likes in 397 posts.

Group: User

well it's just a matter of time that the DEV, with Radovan and others invaluable support and clever betatesting, gets to include great features like nonlinear solvers, symbolic calculations of every kind, new ways of documenting your results and powerful graphing capabilities.
Maybe it would help if someone could elaborate a step by step tutorial on how to create a SMath plugin taking any of the math libraries available on the net
#5 Posted: 5/24/2012 7:53:21 AM
Davide Carpi

Davide Carpi

1415 likes in 2872 posts.

Group: Moderator

Hi,

i've an issue with solve() and functions...

This is a part of a script to calculate properties for generic polygons (PolyProperties_DEBUG.zip

Looking for the axes that split in two chunks with equal area the polygon, instead of an iterative procedure i've thought to use solve(), but doesn't work



anyone have an idea to fix it?


regards,

w3b5urf3r
If you like my plugins please consider to support the program buying a license; for personal contributions to me: paypal.me/dcprojects
#6 Posted: 5/24/2012 3:27:07 PM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Hello w3b5urf3r,

It seems that solve() and roots() will not work with eval() sometimes. Here is a simple example:
solve(), roots(), eval()
You used eval() in two of your functions in the same way like f(x):=eval(...). I tried to remove eval() from both of them but even then solve() did not work. Although you can get single values by calculating root function, there can not be made a plot of the functions you are searching roots from - with or without eval. I do not know why. Any explanation that comes to my mind is that there may be some internal problem of the solve() function in your example. Maybe some kind of bug - not sure about it.

I've noticed another problem in resolving this. As we can calculate values from the root function, I get few values and made independent and dependent vectors. After that, I used interpolation and this way solve() will work. The strange thing is that there are quite different results with and without eval() used. Here is the part of your corrected file - using eval:
witheval
PolyProperties_DEBUG-witheval.sm

And without eval()
withouteval
PolyProperties_DEBUG-withouteval.sm

As you can see, there is quite a difference. The only difference between these two files is in using eval(). This might be a bug I have no idea where it comes from. I suppose that the version with eval() is the correct one.

As a final solution to your problem, if these functions are linear then you do not need solve(). Just calculate the intersection with x axis from a straight line (two points would be enough for that).

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#7 Posted: 5/28/2012 7:00:45 AM
Davide Carpi

Davide Carpi

1415 likes in 2872 posts.

Group: Moderator

Hi omorr

Yes I agree with you, there's something strange...

the second solution is the correct one; unfortunately for my purposes the functions are not linear, so I've built a dedicated script (I've tried with the secant method but sometimes diverges)

I'll publish the entire script just finished debugging


regards,

w3b5urf3r
If you like my plugins please consider to support the program buying a license; for personal contributions to me: paypal.me/dcprojects
#8 Posted: 5/28/2012 10:32:59 AM
Edward Ulle

Edward Ulle

20 likes in 182 posts.

Group: Moderator

Radovan,

On your original topic of algorithms for SMath, its not as easy as plugging in a module from another author. SMath does not store numbers as floating point but a series of Terms. For example 2/3 is not stored as 0.6666...7 but 2 and 3 and /. I believe Andrey did it this way for Optimization and it turns out worked well for storing units. To really take advantage of the SMath engine an algorithm would have to be reprogrammed to use the SMath object library.

If I had more time and was a better programmer I'd try to take on some of these. But alas I'm a trial and error programmer. Lots of errors but eventually get it right. Maybe a suggestion for an easy one and I'll try.
Ed
1 users liked this post
Radovan Omorjan 5/28/2012 11:14:00 AM
#9 Posted: 5/28/2012 11:23:16 AM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Thank you Ed for your reply and your good will

Regarding your answer I am a bit, let me say, discouraged. I just thought (maybe naively) that some numerical library could be incorporated into SMath (with more or less programming effort) and that someone will find the time and make that sooner or later. If I understood you well, that would not be an easy task - maybe impossible. Never mind, I just hope that some kind of quite powerful nonlinear equation solvers will be available for SMath sooner or later, among other numerical stuff. At the moment solve() and roots() will do, but unfortunately they are not good enough.

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#10 Posted: 6/3/2012 2:19:42 PM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Hello,

Here is a simple realization of Broydes's method (Quasi-Newton) for solving the system of nonlinear equations.

Broydens' method

Broyden-method-2.sm

Try it when rots() fail. Moreover, there will be no complain if the functions were defined with eval().

Regards,
Radovan

When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
1 users liked this post
Andrey Ivashov 6/3/2012 4:56:00 PM
#11 Posted: 6/4/2012 3:17:41 PM
kilele

kilele

133 likes in 397 posts.

Group: User

Is this infallible?
I suppose that SMath's solve() investigates the particular function entered by the user and try different strategies/algorithms/intervals.
If I remember, the Bisection method always converge thanks to Bolzano theorem, i.e provided that the function is continuous and there's a change of sign between the two points. Bisection is slow but who cares (Soyuz spacecrafts are ugly but also reliable) if it manages to find a bunch of roots combined with other algorithms and strategies to seek roots within certain intervals based on the kind of function.
#12 Posted: 6/4/2012 5:26:46 PM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Hello,
Wrote

Is this infallible?

No, it is not of course. It is just an alternative to roots(). It worked for me for some problems when roots() failed

Wrote

If I remember, the Bisection method always converge ...

You are right. I do not know how solve() works. But Bisection method for n-dimensional problem in SMath? Hmm...I am not so sure about that.

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#13 Posted: 6/4/2012 8:54:29 PM
kilele

kilele

133 likes in 397 posts.

Group: User

Wrote

But Bisection method for n-dimensional problem in SMath? Hmm...I am not so sure about that.



There seems to be bibliography about generalized bisection, of course this is advanced maths far beyond my grasp :d
I might be able to access some of these articles, tell me by pm if someone is interested.

http://www.ddm.org/DD07/Mejzlik.pdf
http://www.springerlink.com/content/w72615872r512112/
http://dl.acm.org/citation.cfm?id=2705
http://dl.acm.org/citation.cfm?id=11151
http://goo.gl/3uKT4
http://goo.gl/qzPPB
http://goo.gl/2D4E2
1 users liked this post
Radovan Omorjan 6/5/2012 2:25:00 AM
#14 Posted: 6/5/2012 4:12:40 AM
Radovan Omorjan

Radovan Omorjan

325 likes in 2052 posts.

Group: Moderator

Thanks kilele for the reply,

Actually, we could try to make all those functions in SMath, but I think there is no point to make them ourselves (except for our personal use), like this one:

bisect
bisect.sm

By the way, this will not work if we use eval(), as well.

eval()

Do not know why is this not working

We could also try to make an n-dimensional bisect as well like an SMath function (with some more effort). But this is not the point. I think SMath should have those functions already incorporated, IMHO. The problem is how to make them. I thought the more practical way is to use some already available code and make some plugins from them - but it seems I was wrong.

Regards,
Radovan
When Sisyphus climbed to the top of a hill, they said: "Wrong boulder!"
#15 Posted: 6/5/2012 5:54:15 AM
kilele

kilele

133 likes in 397 posts.

Group: User

Thanks again for your valuable contributions..

I've found an article at
http://numericalmethods.eng.usf.edu
explaining the bisection method and its advantages/drawbacks with examples:
http://numericalmethods.eng.usf.edu/mws/gen/03nle/mws_gen_nle_txt_bisection.doc
#16 Posted: 9/5/2012 7:46:25 AM
kilele

kilele

133 likes in 397 posts.

Group: User

a question for the experts? what is "homotopy method" ?
It seems to be a very successful method in solving the roots for polynomial systems. Look at the end of this document
  • New Posts New Posts
  • No New Posts No New Posts