Moving sofa demo script

[GoogleCodeExporter]

It would be awesome with a demo script that calculates the constant in
Gerver's solution to the moving sofa problem with high precision:

http://mathworld.wolfram.com/MovingSofaProblem.html

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.

Original issue reported on code.google.com by fredrik....@gmail.com on 27 Aug 2008 at 9:48

[GoogleCodeExporter]

Indeed. These movers problems are always interesting and mostly nontrivial.

Original comment by ondrej.c...@gmail.com on 28 Aug 2008 at 12:03

[GoogleCodeExporter]

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.

Original comment by Vinzent.Steinberg@gmail.com on 28 Aug 2008 at 10:06

[GoogleCodeExporter]

The result is not exactly the same, I suspect there is a typo or bug somewhere. 

Anyway, have fun!

For the impatient:

$ python sofa.py
[...]
step 5:
x0: [ 0.109131677465481]
[  1.39296704533563]
[0.0453589103507513]
[  1.02936055239276]
fx: [1.33573707650214e-16]
[5.14779191496118e-16]
[                   0]
[1.66533453693773e-16]
||fx||: 5.14779191496118e-16
Jx:
[-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]
[2.31251759221299e-16]
[2.83086607741439e-17]
[1.02103115770459e-16]
canceled, won't get more excact
[ 0.109131677465481]
[  1.39296704533563]
[0.0453589103507513]
[  1.02936055239276]

Original comment by Vinzent.Steinberg@gmail.com on 28 Aug 2008 at 10:36

Attachments:

• sofa.py

[GoogleCodeExporter]

Got this when trying to start from (0., 0., 0., 0.):

Traceback (most recent call last):
  File "sofa.py", 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/solvers.py", line 602,
in msolve
    x = newton(f, x0, J, **kwargs)
  File "/usr/local/lib/python2.5/site-packages/sympy/solvers/numeric.py", line 52, in
newton
    fx = f(*x0)
  File "<string>", line 1, in <lambda>
  File "/usr/local/lib/python2.5/site-packages/sympy/mpmath/mptypes.py", line 646, in f
    x = convert_lossless(x)
  File "/usr/local/lib/python2.5/site-packages/sympy/mpmath/mptypes.py", line 195, in
convert_lossless
    raise TypeError("cannot create mpf from " + repr(x))
TypeError: cannot create mpf from 0

Original comment by Vinzent.Steinberg@gmail.com on 28 Aug 2008 at 10:40

[GoogleCodeExporter]

There is a phi instead of a theta in the second equation.

Fixing it gives:

[0.0944265608436529]
[  1.39920372733355]
[0.0391773647900836]
[ 0.681301509382725]

Nice!

But I already knew about sympy.msolve. I would like to have a function like 
this in
mpmath, not requiring any symbolic processing.

Original comment by fredrik....@gmail.com on 28 Aug 2008 at 11:15

[GoogleCodeExporter]

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 numeric.py in 
the
solvers directory. It can easily be moved to mpmath, only msolve has to be 
fixed to
find it.

Original comment by Vinzent.Steinberg@gmail.com on 29 Aug 2008 at 8:13

[GoogleCodeExporter]

This example should work now in mpmath, using findroot().

Original comment by Vinzent.Steinberg@gmail.com on 16 Jan 2009 at 8:21

• Changed state: Started

[GoogleCodeExporter]

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])
    matrix(
    [['-0.942845232974915'],
     ['0.68171791938912']])
    >>> nsolve([e1,e3],[Phi,theta],[0.2,.5])
    matrix(
    [['0.0453589103507513'],
     ['1.02936055239276']])

Original comment by smi...@gmail.com on 9 Dec 2011 at 4:55