# endless loop using ode_solve

 0 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? asked Aug 22 '10 ngativ 93 ● 1 ● 7 ● 15 Just to check: Are you solving Legendre's DE? Mitesh Patel (Aug 22 '10) 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 (Aug 22 '10) how i do backtrace? the script doesn't crash.... it just calculate forever ngativ (Aug 23 '10) 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 (Aug 23 '10)

 0 It seems like this is just ridiculously slow, but I'm not sure why. I tried tweaking your code in a variety of ways and nothing helped. To try to figure out what was going on I wrote the following Runge-Kutta 4th order code in Cython; its pretty hackish but it makes it clear that the solution isn't doing anything very strange. The version below is hard-coded for your ll=2. %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  and then you can get a plot with: ans = RK4_2d(leg_f1, leg_f2, 0.0, 0.0, 1.0, 1.0, 40) list_plot([[x[0],x[1]] for x in ans], plotjoined = True)  Anyway for me it seems easier to implement such things directly than to wrestle with the GSL wrapping. posted Aug 27 '10 mhampton 359 ● 2 ● 6 ● 15 is a convergence problem for l!=1. ngativ (Aug 28 '10)

[hide preview]

## Stats:

Asked: Aug 22 '10

Seen: 123 times

Last updated: Aug 27 '10