endless loop using ode_solve

asked 2010-08-22 23:26:49 +0200

ngativ
103 ●6 ●11 ●16

updated 2010-08-23 00:55:41 +0200

Im solving a system equation within a for loop, but it never ends But its ok to evaluate f1(n) outside the loop.... why?

def f1(ll):

    def f1(ll):
        T=ode_solver()
        f = lambda t,y:[y[1],eqa.subs(l=ll,s=y[0],w=y[1],x=t)]
        T.function=f
        T.y_0=[0,1]
        T.algorithm="rk4"
        T.t_span=[0,1]
        T.h=1
        T.ode_solve(num_points=40)
        s=T.solution[-1][1][0]-1
        return s

for i in range(3) :f1(i)


doesn't works except when f(1) ,otherwise  the solver is unable to end. But f(1) or eqa(l=1) is the real solution--- :S

where eqa is

var('l s x w')
eqa=((l^2 + l)*s - 2*w*x)/(x^2 - 1)

in the other hand this works

def g(ll):
    Sol=desolve_system_rk4([w,eqa(l=ll)],[s,w],ics=[0,0,1],end_points=[1.0],ivar=x,step=0.01)
    return (Sol[len(Sol)-1][1]-1)
for i in range(3) :f1(i)

so.. whats the problem with ode_solver?

edit retag flag offensive close merge delete

Comments

Just to check: Are you solving Legendre's DE?

Mitesh Patel ( 2010-08-23 01:50:47 +0200 )edit

yea, im just trying to learn Sage solving problems like this. for Pl(x)=x i found the eigenvalue l= 0.994838832387 using desolve_system_rk4. I just wanna to solve the problem using ode_solver and odeint.

ngativ ( 2010-08-23 01:57:29 +0200 )edit

how i do backtrace? the script doesn't crash.... it just calculate forever

ngativ ( 2010-08-23 02:00:59 +0200 )edit

Well... it seems to be that when the solver evaluates the function (eqa(t=1)) the solution diverges because the 1/t factor. But this doesn't happens with the desolver_system_rk4. i dont like this, maybe ode_solver is faster but.. it really has to have this kind of problems?

ngativ ( 2010-08-23 05:31:13 +0200 )edit

add a comment

%cython cdef double ll = 2.0 cpdef double leg_f1(double t, double y, double yp): return yp cpdef double leg_f2(double t, double y, double yp): return ((ll**2 + ll)*y - 2.0*yp*t)/(t**2 - 1.0) cpdef RK4_2d(f1, f2, double t_start, double y_start, double yp_start, double t_end, int steps): cdef double step_size = (t_end - t_start)/steps cdef double t_current = t_start cdef double y_current = y_start cdef double yp_current = yp_start cdef list answer_table = [] cdef int j answer_table.append([t_current,y_current, yp_current]) for j in range(0,steps): k1=f1(t_current, y_current, yp_current) k2=f1(t_current+step_size/2, y_current + k1*step_size/2, yp_current + k1*step_size/2) k3=f1(t_current+step_size/2, y_current + k2*step_size/2, yp_current + k1*step_size/2) k4=f1(t_current+step_size, y_current + k3*step_size, yp_current + k1*step_size/2) y_current = y_current + (step_size/6)*(k1+2*k2+2*k3+k4) k1=f2(t_current, y_current, yp_current) k2=f2(t_current+step_size/2, y_current + k1*step_size/2, yp_current + k1*step_size/2) k3=f2(t_current+step_size/2, y_current + k2*step_size/2, yp_current + k1*step_size/2) k4=f2(t_current+step_size, y_current + k3*step_size, yp_current + k1*step_size/2) yp_current = yp_current + (step_size/6)*(k1+2*k2+2*k3+k4) t_current += step_size answer_table.append([t_current,y_current,yp_current]) return answer_table

Comments

is a convergence problem for l!=1.

ngativ ( 2010-08-28 13:31:36 +0200 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

endless loop using ode_solve

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

endless loop using ode_solve edit

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

endless loop using ode_solve