I think you have the function slightly wrong after looking at the
paper. I tried the following - the cython code has to be in a
separate cell from the second part:
%cython
cdef double f(double z0, double z1, double x, double b):
return (z1**2+1)*z1/x + (2/b-z0)*(1+z1**2)**(1.5)
cpdef RK4_drop(double t_start, double z0_start, double z1_start,
double t_end, int steps, double b):
'''
Fourth-order Runge-Kutta solver for drop model.
'''
cdef double step_size = (t_end - t_start)/steps
cdef double t_current = t_start
cdef double z0_current = z0_start
cdef double z1_current = z1_start
cdef double k10,k20,k30,k40,k11,k21,k31,k41
cdef list answer_table = []
cdef int j
answer_table.append([t_current,z0_current,z1_current])
for j in range(0,steps):
k10 = z1_current
k11 = f(z0_current, z1_current, t_current, b)
k20 = z1_current + k11*step_size/2
k21 = f(z0_current + k10*step_size/2, z1_current +
k11*step_size/2, t_current, b)
k30 = z1_current + k21*step_size/2
k31 = f(z0_current + k20*step_size/2, z1_current +
k21*step_size/2, t_current, b)
k40 = z1_current + k31*step_size
k41 = f(z0_current + k30*step_size, z1_current +
k31*step_size, t_current, b)
t_current += step_size
z0_current = z0_current + (step_size/6)*(k10+2*k20+2*k30+k40)
z1_current = z1_current + (step_size/6)*(k11+2*k21+2*k31+k41)
answer_table.append([t_current,z0_current,z1_current])
return answer_table
and then I tried two values for B:
traj18 = RK4_drop(.1,.1,0.01,.5,300,1.8)
traj20 = RK4_drop(.1,.1,0.01,.5,300,2.0)
list_plot([x[0:2] for x in traj18])+list_plot([x[0:2] for x in
traj20])
I don't understand the model well enough to tell if this is correct.
Hope that helps,
M. Hampton
On Dec 2, 11:11 am, ynk1 <[email protected]> wrote:
> 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
