My favorites | Sign in
Project Home Downloads Wiki Issues Source
New issue   Search
  Advanced search   Search tips   Subscriptions
Issue 56: Moving sofa demo script
3 people starred this issue and may be notified of changes. Back to list
Status:  Started
Owner:  ----

Sign in to add a comment
Project Member Reported by, Aug 27, 2008
It would be awesome with a demo script that calculates the constant in
Gerver's solution to the moving sofa problem with high precision:

The main missing component is a solver for the 4x4 nonlinear problem
(Vinzent, how far off is this?).

Then there are just some double integrals that should be possible to
calculate already, although adaptive quadrature might be necessary due to
their piecewise nature.
Aug 28, 2008
Project Member #1
Indeed. These movers problems are always interesting and mostly nontrivial.
Aug 28, 2008
Project Member #2
The solver exists and works, your starting point just has to be close enough to the
root. Depending on the system, "close" can mean +-0.001. It converges only locally,
it needs a special strategy to converge globally. Sadly this is much more complicated
than the current algorithm and beyond my math skills.
Aug 28, 2008
Project Member #3
The result is not exactly the same, I suspect there is a typo or bug somewhere. 

Anyway, have fun!

For the impatient:

$ python
step 5:
x0: [ 0.109131677465481]
[  1.39296704533563]
[  1.02936055239276]
fx: [1.33573707650214e-16]
[                   0]
||fx||: 5.14779191496118e-16
[-0.483604523522736, -0.0906867162958127,  -2.33985908666742, -0.0798124248931422]
[  2.61625220092966,   -1.99794292197937, -0.992947564121786,   -2.26615019314999]
[ 0.998971460989687, -0.0453433581479064,  -2.41812586118966,                   0]
[    1.492000821021,                  -1,  -2.04656665975375,  0.0465666597537454]
s: [7.90209121450792e-17]
canceled, won't get more excact
[ 0.109131677465481]
[  1.39296704533563]
[  1.02936055239276]
599 bytes   View   Download
Aug 28, 2008
Project Member #4
Got this when trying to start from (0., 0., 0., 0.):

Traceback (most recent call last):
  File "", line 17, in <module>
    print msolve((A, B, phi, theta), system, (0., 0., 0., 0.), verbose=1)
  File "/usr/local/lib/python2.5/site-packages/sympy/solvers/", line 602,
in msolve
    x = newton(f, x0, J, **kwargs)
  File "/usr/local/lib/python2.5/site-packages/sympy/solvers/", line 52, in
    fx = f(*x0)
  File "<string>", line 1, in <lambda>
  File "/usr/local/lib/python2.5/site-packages/sympy/mpmath/", line 646, in f
    x = convert_lossless(x)
  File "/usr/local/lib/python2.5/site-packages/sympy/mpmath/", line 195, in
    raise TypeError("cannot create mpf from " + repr(x))
TypeError: cannot create mpf from 0

Aug 28, 2008
Project Member #5
There is a phi instead of a theta in the second equation.

Fixing it gives:

[  1.39920372733355]
[ 0.681301509382725]


But I already knew about sympy.msolve. I would like to have a function like this in
mpmath, not requiring any symbolic processing.
Aug 29, 2008
Project Member #6
In fact, it's based on a numerical function. The only symbolic math is done to
calculate the Jacobian, and this can be done numerically (it's even faster for really
complicated functions). I just did not yet implement it, but it's trivial.

I separated the numerical stuff, you might want to have a look at in the
solvers directory. It can easily be moved to mpmath, only msolve has to be fixed to
find it.
Jan 16, 2009
Project Member #7
This example should work now in mpmath, using findroot().
Status: Started
Dec 9, 2011
Just a reminder that the system can be reduced to two equations if you eliminate A and B from the system. (e.g. solve expressions 0 and 2 for A and B and then solve the remaining two equations for Phi and theta). Using sympy's nsolve gives:

    >>> nsolve([e1,e3],[Phi,theta],[0.1,.5])
    >>> nsolve([e1,e3],[Phi,theta],[0.2,.5])

Sign in to add a comment

Powered by Google Project Hosting