My favorites | Sign in
Project Home Downloads Wiki Issues Source
New issue   Search
for
  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 fredrik....@gmail.com, 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:

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.
Aug 28, 2008
Project Member #1 ondrej.c...@gmail.com
Indeed. These movers problems are always interesting and mostly nontrivial.
Aug 28, 2008
Project Member #2 Vinzent.Steinberg@gmail.com
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 Vinzent.Steinberg@gmail.com
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]


sofa.py
599 bytes   View   Download
Aug 28, 2008
Project Member #4 Vinzent.Steinberg@gmail.com
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

Aug 28, 2008
Project Member #5 fredrik....@gmail.com
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.
Aug 29, 2008
Project Member #6 Vinzent.Steinberg@gmail.com
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.
Jan 16, 2009
Project Member #7 Vinzent.Steinberg@gmail.com
This example should work now in mpmath, using findroot().
Status: Started
Dec 9, 2011
#8 smi...@gmail.com
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']])

Sign in to add a comment

Powered by Google Project Hosting