0

I'm trying to solve (for m_0) numerically the following ordinary differential equation:

dm0/dx=(((1-x)*(x*(2-x))**(1.5))/(k+x)**2)*(((x*(2-x))/3.0)*(dw/dx)**2 + ((8*(k+1))/(3*(k+x)))*w**2)

The values of w and dw/dx have been found already numerically using the Runge-Kutta 4th order and k is a factor that is fixed. I wrote a code where I call the values for w and dw/dx from an external file, then I organize them in an array, then I call the array in the function and then I run the integration. My outcome is not what it's expected :(, I don't know what is wrong. If anyone could give me a hand, it would be highly appreciated. Thank you!

from math import sqrt
from numpy import array,zeros,loadtxt
from printSoln import *
from run_kut4 import *

m = 1.15                         # Just a constant.
k = 3.0*sqrt(1.0-(1.0/m))-1.0    # k in terms of m.

omegas = loadtxt("omega.txt",float)    # Import values of w
domegas = loadtxt("domega.txt",float)  # Import values of dw/dx

w = []    # Defines the array w to store the values w^2
s = 0.0
for s in omegas:
    w.append(s**2)       # Calculates the values w**2
omeg = array(w,float)    # Array to store the value of w**2

dw = []      # Defines the array dw to store the values dw**2
t = 0.0
for t in domegas:
    dw.append(t**2)    # Calculates the values for dw**2
domeg = array(dw,float)     # Array to store the values of dw**2

x = 1.0e-12              # Starting point of integration
xStop = (2.0 - k)/3.0    # Final point of integration 

def F(x,y):      # Define function to be integrated
    F = zeros(1)
    for i in domeg:  # Loop to call w^2, (dw/dx)^2       
        for j in omeg:
            F[0] = (((1.0-x)*(x*(2.0-x))**(1.5))/(k+x)**2)*((1.0/3.0)*x* (2.0-x)*domeg[i] + (8.0*(k+1.0)*omeg[j])/(3.0*(k+x))) 
            return F

y = array([((32.0*sqrt(2.0)*(k+1.0)*(x**2.5))/(15.0*(k**3)))])  # Initial condition for m_{0}

h = 1.0e-5          # Integration step
freq = 0            # Prints only initial and final values
X,Y = integrate(F,x,y,xStop,h)      # Calls Runge-Kutta 4
printSoln(X,Y,freq)                 # Prints solution
CamPos
  • 1
  • Please go back to basic syntax and semantics. `for i in domeg` iterates over the elements of `domeg`, `domeg[i]` is thus non-sensical. The return statement is misplaced, the double loop in itself is non-sensical, since `F[0]` only contains the result of the last performed iteration. ... Then we can start on how to solve coupled ODE systems. – Lutz Lehmann Sep 08 '15 at 13:39
  • I imagined my loops are wrong, it was just a trial and error thing. I'm a beginner in Python and I'm sort of in a rush to solve this. Thanks. – CamPos Sep 08 '15 at 14:55

1 Answers1

0

Interpreting your verbal description, there is an ODE for omega, w'=F(x,w), and a coupled ODE for m0, m'=G(x,m,w,w'). The almost always optimal way to solve this is to treat it as system of ODE,

def ODEfunc(x,y)
    w,m = y
    dw = F(x,w)
    dm = G(x,m,w,dw)
    return np.array([dw, dm])

which you can then insert in the ODE solver of your choice, e.g., the fictitious 

ODEintegrate(ODEfunc, xsamples, y0)
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Hey, thanks for the reply. The ODE for w was already solved numerically, and these values are correct. I got values for w and (dw/dx) in the same range as above, using the same step. As you can see, this equation for m_0 involves (dw/dx)^2 and w^2. Clearly if you have some function w(x), is easy, you just put into the equation, take the square, integrate and done. But here, I have numbers for w and (dw/dx) so I don't understand how to incorporate them here to solve the ODE for m0. That's what I was trying to do in that incorrect attempt with the double for loop. – CamPos Sep 08 '15 at 15:05
  • Assuming that the step sizes of the integration of w and m0 are the same, you will at least have to interpolate the half-step values required in the RK4 scheme. Then again it is a design question of how to transmit the index into the w arrays or better the w and dw/dx values to the equation for m0. – Lutz Lehmann Sep 08 '15 at 16:43