hi,
I am a sage newbie willing to solve numerically a non-linear ode: the
equation for a hanging liquid drop from a capillary.
it reads
d2z/dt2 1 dz/dt
------------------------------ + ---- -----------------------------
+ z - 2/b = 0
( 1 + (dz/dt)^2 )^(3/2) x ( 1 + (dz/dt)^2 )^(1/2)
where z in a function of x and b is a parameter.
The problem is describe in the paper: Meister & Latychevskaia,
Journal of Chemical Education, vol83 p117 (2006).
The apex of the drop is at x=0 z(0)=0 and symmetric with the vertical
z-axis
They integrate it with a "4th order Runge-Kutta interpolation".
Following various examples, isolating d2z/dt2, I have tried naively
this code in the notebook:
import scipy
from scipy import integrate
def func(z,x,b):
return [ z[1], (z[1]**2+1)*z[1]/x + (2/1-
z[0]-1)*(1+z[1]**2)**(3/2) ]
def jacf(y,t):
return [ [0,1], [-(z[1]**2 + 1.)**(3/2), -3*sqrt(z[1]**2 +
1.)*(z[0] - 2/1.)*z[1] + (3.*z[1]**2 + 1.)/x ]]
x = scipy.arange(0, 2, 0.01)
drop = integrate.odeint(func, [0,0], x, args=(1,), Dfun=jacf)
which gives me a list of NaN.. I guess because of the 1/x
Changing the range of x for curiosity, the system crashes... in the
sage terminal I get:
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] got EOF subprocess
must have crashed...
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] At line 84 of file
scipy/integrate/odepack/xerrwv.f (Unit 6)
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] Traceback: not
available, compile with -ftrace=frame or -ftrace=full
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] Fortran runtime
error: Expected CHARACTER for item 1 in formatted transfer, got
INTEGER. If you want to make character descriptors typeless, compile
with -fsloppy-char
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] (1x,15a4)
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1] ^
2010-12-02 16:31:53+0100 [HTTPChannel,0,127.0.0.1]
and I have to restart the worksheet. I may even have to quit and come
back again otherwise I get various NameError message like "scipy not
defined" "drop not defined", when evaluating cells and even loose the
<tab> completion. Changing some integers for floats in func or drop is
not changing anything.
My poor knoweldge of sage/python, and math in general, is certainly
not helping me, So..
- What/where is the problem ?
- is there a better way than with scipy's odeint?
- How to define the problem to deal with x=0? if that is the
problem...
- or do I get it totally wrong... :-\
- Also for small b values the drop shape should have a neck... meaning
that several z exist for one value of x and there is x0 such that (dz/
dx)_x0 = +oo
it seems that the definition of x in the code is not appropriate for
that...
Any help on sage and/or the math would be very much appreciated.
Yves
sage 4.3 + opensuse 11.1
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
