11

I posted some benchmarks here: http://www.cecm.sfu.ca/~rpearcea/mgb.html (archived copy) These are for total degree orders. To solve systems you typically need to do more work. Timings are for a typical midrange desktop as of 2015 (Haswell Core i5 quad core). The fastest system on one core is Magma, which uses floating point arithmetic and SSE/AVX. Magma ...


10

I don't think I can answer all of these, but I can give you some thoughts from what I have seen. I think symbolic computation is a great tool, but it obviously most useful for problems that have analytical solutions and forms. Not every problem is really represented this way, however. So this is where numeric computation would probably prove a better option....


9

It's state-dependent control flow that's an issue. function f(u) z = 0.0 while z < 10.0 z += u u += z^2 end return u end What's the program for computing the derivative? Automatic differentiation would give you: function f(_u) u = (_u,1.0) # seed the input derivative for the jvp in direction of basis e1 z = (0.0,0.0) while z[1] &...


8

You can solve this numerically in Python without symbolic computation. from __future__ import print_function, division import numpy as np from numpy import exp from scipy.integrate import quad from scipy.optimize import root def f1(a1, a2, x): return exp(a1 * x + a2 * x * x * x) / (1 + x * x) def f2(a1, a2, x): return exp(a1 * x + a2 * x * x * x) *...


8

Googling benchmark polynomial systems leads to a few hits, including the University of Mannheim's Computer Algebra Benchmark Initiative. Sadly, most of these are out of date or defunct. The most active seems to be the SymbolicData Wiki, but as far as I can tell, it only collects benchmark problems, not benchmark results. Some comparisons (dating back to ...


6

Always keep in mind that the results of any benchmark will depend, in addition to the problem's size, on the base field over which the polynomial ring is defined (rational numbers or integers modulo some power of a prime number). The FGb library is an actively developed and high-performance implementation of the F5 algorithm. A benchmark comparing FGb to ...


6

Before answering much, I want to point out something very important that this question sounds like it should touch on: the fact that numeric computation is kind of meaningless without a symbolic part. You don't just dump your problem onto a numeric solver and let it chug unless it's a very, very simple one (even then, someone spent a lot of time at the ...


6

One word: Modularity. There are a lot of repeated expressions in your Jacobian that could be written as their own function. There's no reason for you to write the same operation more than once, and that will make debugging easier; if you only write it once there's only one place for an error (in theory). Modular code will also make testing easier; ...


6

I wrote a Python package called PyPGE. PyPGE is a Symbolic Regression implementation based on Prioritized Grammar Enumeration (1), not Evolutionary or Genetic Programming. It produces a deterministic Symbolic Regression algorithm. (1) Worm, Tony, and Kenneth Chiu. "Prioritized grammar enumeration: symbolic regression by dynamic programming." Proceedings of ...


6

Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale in the dynamics ... You need to use a reversible ODE solver method if you want to do this. I actually recently showed in a blog post that there are many cases ...


5

There are several strategies to consider: Find the derivatives in symbolic form using a CAS, then export code for computing the derivatives. Use an automatic differentiation (AD) tool to produce code that computes the derivatives from code to compute the functions. Use finite difference approximations to approximate the Jacobian. Automatic ...


5

But how do you numerically (Runge-Kutta or whatever) solve the objects coming out of SymPy? That depends a lot on your specific problem and what level of optimization you need. Most ODE solvers (such as those included in SciPy) require that you provide them a Python function representing the right-hand side of the ODE. As you do not want things to be ...


4

Philips in 1988 proved the following relationship: If $f(x)$ is an infinitely differentiable function defined on the interval $[-1,+1]$ and its Legendre expansion is given by $f(x) = \sum_{n=0}^{\infty}a_{n}P_{n}(x)$, then the Legendre coefficients $a_{n}^{(q)}$ of the $q$-th derivative of $f(x)$ are given by $$a_{n}^{(q)} = \frac{(2n+1)}{2^{q-2}(q-1)!} \...


4

I once started writing anopen source version of Eureqa in Java. The project has limited capabilities but it implements the fitness function described in [1] and couple optimizations mentioned by the authors in other publications (e.g., searching for solutions in Pareto front). Link: https://github.com/pkoperek/hubert [1] Schmidt, Michael, and Hod Lipson. "...


4

After a cursory google search on the subject, it appears that "symbolic regression" is a problem that lends itself greatly to stochastic optimization algorithms like genetic programming (GP). It is conceivable that you should look for an open source genetic programming library with modules specifically for symbolic regression, such as DEAP (Distributed ...


4

Using the function matlabFunction increases the efficiency by a lot! syms x f = sin(x); fp = diff(f,x); F = matlabFunction(fp); tic for i=1:1000 F(0); end toc tic for i=1:1000 cos(0); end toc Elapsed time is 0.037475 seconds. Elapsed time is 0.004203 seconds. Furthermore if the result of the symbolic calculations needs to be saved as a ...


4

Sagemath looks fast enough. Here is an example created by taking a 999x1000 random rational matrix and appending its column sums at the bottom, so that [1, 1, ... , 1, -1] is the kernel. sage: A = random_matrix(QQ, 999, 1000) sage: B = A.stack(sum(A)) sage: %time v = B.kernel() CPU times: user 5.21 s, sys: 36.4 ms, total: 5.24 s Wall time: 5.24 s Vector ...


4

The paper you linked answers the question. Autodiff (or hand differentiation) can differentiate branched program statements. For example, limiters, entropy fixes branching in flux statements, and the like. It can be rather helpful for min max statements as well. You can see an example below: Function(Vn_bar, a_bar, ul, cl, ur, cr) lambda1 = abs(...


4

It's hard to get around the allocations implicit to SymPy in this case. It wants to allocate the matrix, so the easiest thing to do would be, as you show, build individual scalar functions. But then composing those together can be a bit of a hassle, since you don't want to put them into an array since they are all different types and that would then ruin the ...


3

I found the gramEvol R package flexible and easy to use. They have a small tutorial in which they rederive Kepler's third law from data. Note that it relies on Genetic Programmic for its optimisation and thus might return different results if you run it twice.


3

Here is an example of where we have used automatic differentiation using Sacado in one code: http://www.dealii.org/developer/doxygen/deal.II/step_33.html


3

Having used Mathematica, then trying Sage, and now SymPy/SymEngine in Julia, the clear winner is SymEngine.jl. Sage was hard to get working, was very slow, and I found it very hard to develop my own algorithms (which would get decent performance. That's the key: you can write algorithms for it, but they won't necessary be speedy!). The syntax can be really ...


3

Sage would be perfect for classroom use. I once assisted a professor with classroom demonstrations of Physics, and it worked out pretty well. Working with Sage does not require advanced programming knowledge. Students would be able to do it without breaking their heads like in a usual programming class. You can get started with this tutorial or this book. ...


3

tl;dr: Use MATLAB primarily for numerical computation, not symbolic computation. MATLAB was primarily designed as a numerical computation package, and its symbolic capabilities were added later. Also, as a general rule, it's easier to compute quantities with concrete values than it is to calculate the same quantities in the abstract. For instance, no closed ...


3

While I do not know of a toolbox which fits the bill, a good open source alternative to MATLAB and Octave which does have a good solution to this problem is Julia. The linear algebra syntax from MATLAB/Octave almost transfers over to Julia directly, though you need to swap indexing like A(i) to A[i]. But after a quick translation, you can use Julia's ...


3

If you want to find the equilibrium points you need to write your system as a first-order one $$\frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x})$$ and solve the non-linear system of equations $$\mathbf{F}(\mathbf{x})=\mathbf{0}\ ,$$ that might not be solvable analytically. After that, you need to find the eigenvalues of the Jacobian matrix around each ...


3

I want to point out that PyODESys is a Python library that does exactly what you're looking for: send SymPy expressions to ODE solvers. It's nice because it links to more solvers, but @Wrzlprmft's JITCODE is probably better because it compiles the resulting expression and the runtime of the user's function is very important for solving an ODE efficiently.


3

You might look up "complete fraction-free factorization" methods. The paper "Generalized fraction-free LU factorization for singular systems with kernel extraction" contains pseudo-code.


2

A working link is provided in this site: http://vaopt.math.uni-goettingen.de/software.php


Only top voted, non community-wiki answers of a minimum length are eligible