runge kutta 2 in python

Question

I am trying to solve an equation in fluid mechanics using the runge-kutta 2 method, usually it seems quite doable but in this case its with x y and z and i cant seem to make the code. Here is what i have done so far: It does return a plot but I cannot say whether it is at all meaningfull

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Define function to compute velocity field

A = 0.2
B = 1.0
C = 0.20


def velfun(x,y,z,t):
    
    xvelocity =  (A*np.sin(z) +  C*np.cos(x))
    yvelocity =  (B*np.sin(x) + A*np.cos(z))
    zvelocity =  (C*np.sin(y) + B*np.cos(x))
    return xvelocity, yvelocity, zvelocity


# Set stopping time
# and the step length
    
Nstep = 10000   
h     = 0.01

# Create arrays to store x and t values

x     = np.zeros(Nstep);
y     = np.zeros(Nstep);
z     = np.zeros(Nstep);

# Set the initial condition

x[0] = 1,
y[0] = 0,
z[0] = 0

# Carry out steps in a loop

for k in range(1,Nstep):
    
    # Provisional Euler step
    
    tt    = t[k-1]
    xx    = x[k-1]
    yy    = y[k-1]
    zz    = z[k-1]
    ux,uy,uz = velfun(xx,yy,zz,tt)
    xp    = x[k-1] + h*ux
    yp    = y[k-1] + h*uy
    zp    = z[k-1] + h*uz
    # Compute velocity at provisional point
    
    uxp,uyp,uzp = velfun(xp,yp,zp,tt+h)

    # Compute average velocity
    
    uxa = 0.5*(ux + uxp)
    uya = 0.5*(uy + uyp)
    uza = 0.5*(uz + uzp)
    
    # Make final Euler step
    # using average velocity
    
    x[k] = x[k-1] + h*uxa
    y[k] = y[k-1] + h*uya
    z[k] = z[k-1] + h*uza
    


    
# Exact solution
    


# Plot results
fig = plt.figure()
ax = plt.axes(projection='3d')
ax = plt.axes(projection='3d')
ax.scatter3D(x,y,z,'b',label='x (with RK2)')

plt.show()

Answer 1

It is very likely that you want to implement the Arnold-Beltrami-Childress ABC dynamic. Note that the ABC and xyz triples are moved cyclically through each position. This is violated in your code in the last position where you have xzx instead of yzx.

Apart from some small errors or over-corrections (No commas as command separator, as that constructs lists. If the array t is used it must be defined before.) the other parts of your code are correct. Applying all the mentioned points gives a realistic result, extending the time span will make it visibly chaotic.

stereopair of the corrected solution