Issue 56:  Moving sofa demo script  
3 people starred this issue and may be notified of changes.  Back to list 
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
Aug 28, 2008
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
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.33573707650214e16] [5.14779191496118e16] [ 0] [1.66533453693773e16] fx: 5.14779191496118e16 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.90209121450792e17] [2.31251759221299e16] [2.83086607741439e17] [1.02103115770459e16] canceled, won't get more excact [ 0.109131677465481] [ 1.39296704533563] [0.0453589103507513] [ 1.02936055239276]
Aug 28, 2008
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/sitepackages/sympy/solvers/solvers.py", line 602, in msolve x = newton(f, x0, J, **kwargs) File "/usr/local/lib/python2.5/sitepackages/sympy/solvers/numeric.py", line 52, in newton fx = f(*x0) File "<string>", line 1, in <lambda> File "/usr/local/lib/python2.5/sitepackages/sympy/mpmath/mptypes.py", line 646, in f x = convert_lossless(x) File "/usr/local/lib/python2.5/sitepackages/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
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
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
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]) matrix( [['0.942845232974915'], ['0.68171791938912']]) >>> nsolve([e1,e3],[Phi,theta],[0.2,.5]) matrix( [['0.0453589103507513'], ['1.02936055239276']]) 

► Sign in to add a comment 