thanks!!! Helps a lot! Certainly more than I hoped for... Cheers Yves. On Dec 3, 10:42 pm, Marshall Hampton <[email protected]> wrote: > 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
