0

I am somewhat familiar with using RK4 for coupled ODE's, I found a very elegant way (in my opinion) to utilize it in Python, like so:

def RungeKutta4_2(x, y, z, dx, dydx, dzdx):

k1 = dx*dydx(x, y, z)
h1 = dx*dzdx(x, y, z)

k2 = dx*dydx(x+dx/2., y+k1/2., z+h1/2.)
h2 = dx*dzdx(x+dx/2., y+k1/2., z+h1/2.)
k3 = dx*dydx(x+dx/2., y+k2/2., z+h2/2.)
h3 = dx*dzdx(x+dx/2., y+k2/2., z+h2/2.)

k4 = dx*dydx(x+dx, y+k3, z+h3)
h4 = dx*dzdx(x+dx, y+k3, z+h3)

y = y + 1./6.*(k1+2*k2+2*k3+k4)
z = z + 1./6.*(h1+2*h2+2*h3+h4)
x = x + dx

return x, y, z

where x is the independent variable, and y and z are the dependent variables; I like this way of setting it up because the method itself automatically increments x by dx and the function just needs to be called in a while loop; the differential equations themselves are declared in their own functions. For example:

$$ \frac{dy}{dx} = z$$ $$ \frac{dz}{dx} = -y$$

def dydx(x,y,z):
    return z
def dzdx(x,y,z):
    return -y

while x <= x_end:
    x, y, z = RungeKutta4_2(x,y,z,dx,dydx,dzdx)
    x_values.append(x)
    y_values.append(y)
    z_values.append(z)

And then of course y_values and z_values could be plotted against x_values.

Is this the standard way to use RK4 for coupled ODE's? I am currently working on a paper I wish to publish, and my advisor mentioned that it is much more desirable to declare the ODE's in their own array, where dy/dx might be the first element and dz/dx is the second.

I'm having a lot of trouble figuring out why this would be "better", and if "my" way of using RK4 has any noticeable issues that my professor's way would fix?

Vox Winters
  • 1
  • Well, your implementation does seem repetitive - the bulk of it is just doing the same to $z$ (and $h$) what you did to $y$ (and $k$) - so in some sense it's harder to interpret. Further, it is not possible to generalize your method to a system with more variables. The general method I've seen is to take all variables as components of a single vector $X$, then define a single handle $dXdt$ to return the derivatives. That being said, your method looks correct, and I would recommend asking this somewhere else, e.g., StackOverflow. – K. Jiang Feb 25 '24 at 00:16
  • 1
    Thank you for the succinct response! I'll do so, I guess in hindsight this is more of a coding issue than a purely mathematical one. – Vox Winters Feb 25 '24 at 00:36
  • Unless your research is on some ODEs that require writing their own solver, I would recommend using existing implementations. Python isn't great about this, but the solvers in scipy can handle most problems. – whpowell96 Feb 25 '24 at 01:06
  • My needs do require a bit of customization in the solver itself; I figured this was straightforward enough of an implementation that I could apply those customizations easily. – Vox Winters Feb 25 '24 at 01:19
  • You are presenting a borderline case, try the same approach for a system with a state with 10 components and it becomes inconvenient. More on my opinion, and links, in https://math.stackexchange.com/questions/4176684/how-to-implement-a-runge-kutta-method-rk4-for-a-second-order-differential-equa/4177218#4177218 – Lutz Lehmann Feb 25 '24 at 16:21

0 Answers0