coupled ode (URGENT.. PLEASE HELP)

[Samrat Halder]

hi

I am new to python and dont have much experience.

I have a system of coupled first order linear odes:

dy1/dx=f1(x)y1(x) + f2(x)y2(x)

dy2/dx=f2(x)y1(x) +f2(x)y2(x)+(a gaussian function)

y1(1st grid point)=0

y2(last grid point)=0

I have an pre specified x grid and I have the f1, f2, f3, f4 values at all the grid points.

Can I solve this problem using bvp_solver??

I have written the code and I am pasting it over here. But from this one I am not getting the desired nature of y1 and y2.  Please help if there is any logical error in my code. URGENT. Thanks

from __future__ import division

import scikits.bvp_solver

import numpy as np

import numpy as np

import numexpr as ne

from scipy import interpolate

from numpy import array

pi=np.pi

exp=np.exp

x, c, rho, p, Gamma_1, T= np.loadtxt("data1", unpack=True) 

x=x[::-1]

c=c[::-1]

rho=rho[::-1]

p=p[::-1]

Gamma_1=Gamma_1[::-1]

omega=2*pi*0.003

A=1/((2*pi)**0.5)/2.7/10**(-5)

l=500

R=6.955*10**10

sigma2=7.189072*10**(-10)

#defining dx array

x1=np.roll(x,1)

dx=x-x1

dx=dx[1:]

# defining dP array

p1=np.roll(p,1)

dp=p-p1

dp=dp[1:]

# defining drho array

rho1=np.roll(rho,1)

drho=rho-rho1

drho=drho[1:]

x=x[1:]

c=c[1:]

rho=rho[1:]

p=p[1:]

Gamma_1=Gamma_1[1:]

# defining Sl^2 array

#~ Sl2= l*(l+1)*c**2/(R**2)/(x**2)

# defining g array

g=(-1)*dp/dx/rho/R

# N^2 array

N2=(g/R)*(dp/dx/p/Gamma_1 - drho/dx/rho)

#defining the coefficients in the coupled differential equations

c1=R*(g/c**2-2/R/x)

c2=R*(l*(l+1)*c**2/R**2/omega**2/x**2 -1)/(rho*c**2)

d1=R*rho*(omega**2-N2)

d2=R*g/c**2 

#~ z=A*exp((-1)*(((x-1.0+(5.0/7.0)*10**(-4.0))**2)/(2.0*sigma2)))  ##need to be checked

#~ def extrap1d(interpolator):

    #~ xs = interpolator.x

    #~ ys = interpolator.y

#~ 

    #~ def pointwise(x):

        #~ if x < xs[0]:

            #~ return ys[0]+(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])

        #~ elif x > xs[-1]:

            #~ return ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])

        #~ else:

            #~ return interpolator(x)

#~ 

    #~ def ufunclike(xs):

        #~ return np.array(map(pointwise, array(xs)))

#~ 

    #~ return ufunclike

f1=interpolate.interp1d(x,c1)

#f1=extrap1d(f1)

f2=interpolate.interp1d(x,c2)

#f2=extrap1d(f2)

f3=interpolate.interp1d(x,d1)

#f3=extrap1d(f3)

f4=interpolate.interp1d(x,d2)

#f4=extrap1d(f4)

x,t,r,e,w,s=np.loadtxt("data_xgrid", unpack=True)

x=x[::-1]

xileft=0.0

pright=0.0

left=0.5007668 

right=1.0003916

#~ def function(x,y):

#~ 

#~ xiprime=f1(x)*y[0] + f2(x)*y[1]

#~ pprime=f3(x)*y[0] -f4(x)*y[1]+ R*z

#~ 

#~ return np.array([xiprime,pprime])

def function(x,y):

#~ print y

dxi_dx=f1(x)*y[0] + f2(x)*y[1]

dp_dx=f3(x)*y[0] -f4(x)*y[1]+A*R*exp((-1)*(((x-1.0+(5.0/7.0)*10**(-4.0))**2)/(2.0*sigma2)))

return np.array([dxi_dx,dp_dx])

#return yprime

def boundary_conditions(y_left,y_right):

return (np.array([y_left[0]-xileft]),np.array([y_right[1]-pright]))

problem = scikits.bvp_solver.ProblemDefinition(num_ODE = 2,num_parameters = 0,num_left_boundary_conditions = 1,

boundary_points = (left, right),function = function,

boundary_conditions = boundary_conditions)

                                      

                                      

                                    

                                      

                                      

solution = scikits.bvp_solver.solve(problem,solution_guess = ((xileft- pright)/2.0,(xileft- pright)/2.0))

                                             

                            

y= solution(x)

y=y.T

y=np.append(np.atleast_2d(x).T,y,axis=1)

np.savetxt('bvp_solver500', y)

[John Salvatier]

You're not doing anything obviously wrong to me. I would try checking that function returns the values you think it should return.

--

---
You received this message because you are subscribed to the Google Groups "scikits.bvp_solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scikits-bvp_sol...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Samrat Halder]

Hi John
Thanks a lot for the reply. Please look into that and let me know if you found anything wrong. Its related to my summer project and  I really need some help. Thanks again.

You received this message because you are subscribed to a topic in the Google Groups "scikits.bvp_solver" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/scikits-bvp_solver/Hvk5p5QQLUw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to scikits-bvp_sol...@googlegroups.com.