1

I am solving two coupled ordinary differential equations using the 4th-order Runge-Kutta method. I am having trouble in printing the values of z as a result of applying this method. The source code is below for reference. Please, help me fix this code by printing the updated values of z and theta. Thank you.

#Import neeeded modules
import numpy as np
import matplotlib.pyplot as plt

#Input parameters
k = 5 #longitudinal torsional constant
delta = 10**-3 #longitudinal torsional constant
I = 10**-4 #Rotational Inertia
eps = 10**-2 #Coupling constant
m = 0.5


#Time Step
#Setting time array for graph visualization
dt = 0.002 #Time Step
tStop = 0.30 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+dt, dt) #Array of time
z = np.zeros(len(t))
dz = np.zeros(len(t))
theta = np.zeros(len(t))
dtheta = np.zeros(len(t))

#Functions that include the equations of motion
def dYdt(t,u):
    z, dz, theta, dtheta = u
    ddz = -(k*z+0.5*eps*theta)/m
    ddtheta = -(delta*theta+0.5*eps*z)/I
    return np.array([dz, ddz, dtheta, ddtheta])

def rk4(t, u, dt):
    for i in range(len(t)-1):
     # runge_kutta
        k1 = dYdt(t[i], u[i])
        k2 = dYdt(t[i] + dt/2, u[i] + dt/2 * k1)
        k3 = dYdt(t[i] + dt/2, u[i] + dt/2 * k2)
        k4 = dYdt(t[i] + dt, u[i] + dt * k3)
        u.append(u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
        #Unpacking
        z, dz, theta, dtheta = np.asarray(u).T

print(z)

Here are the equations of motions I used to come up with the source code. Coupled ODEs

Rafael Gomez
  • 11
  • 1
  • You are not calling the `rk4` function. I assume you want to do that. And `rk4` doesn't return anything. The last line of `rk4` should be out-dented, so it is outside of the loop. It would be better to return those values, rather than have them be globals, I'll post what I THINK you want below. – Tim Roberts Apr 17 '21 at 03:40
  • @TimRoberts I wanted to call the `rk4` function but I do not know how to do it in that block. Perhaps, this hesitation is reinforced of fear that the code will give so much errors after. Thank you, I will check out your amendments and see if it is what I want. – Rafael Gomez Apr 17 '21 at 03:45
  • Well, MY code doesn't know how to call it, either. It isn't going to run. You need to figure out where t, u, and dt come from. – Tim Roberts Apr 17 '21 at 03:49

1 Answers1

0

This is what I THINK you want, but I don't know what parameters to pass to rk4. Also, u is a 1D vector, so I don't know what you were expecting z, dz, theta, dtheta = np.asarray(u).T to do. As a result, this code won't run, but it will show you a potential design.

import numpy as np
import matplotlib.pyplot as plt

#Input parameters
k = 5 #longitudinal torsional constant
delta = 10**-3 #longitudinal torsional constant
I = 10**-4 #Rotational Inertia
eps = 10**-2 #Coupling constant
m = 0.5


#Time Step
#Setting time array for graph visualization
dt = 0.002 #Time Step
tStop = 0.30 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+dt, dt) #Array of time

#Functions that include the equations of motion
def dYdt(t,u):
    z, dz, theta, dtheta = u
    ddz = -(k*z+0.5*eps*theta)/m
    ddtheta = -(delta*theta+0.5*eps*z)/I
    return np.array([dz, ddz, dtheta, ddtheta])

def rk4(t, u, dt):
    for i in range(len(t)-1):
     # runge_kutta
        k1 = dYdt(t[i], u[i])
        k2 = dYdt(t[i] + dt/2, u[i] + dt/2 * k1)
        k3 = dYdt(t[i] + dt/2, u[i] + dt/2 * k2)
        k4 = dYdt(t[i] + dt, u[i] + dt * k3)
        u.append(u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
    #Unpacking
    return np.asarray(u).T

z, dz, theta, dtheta = rk4(t, u, dt)
print(z)
Tim Roberts
  • 48,973
  • 4
  • 21
  • 30
  • The two coupled odes can be expressed in terms of vector form, which I did in terms of `u`. The code you provided, however, seems to make things a little bit off. – Rafael Gomez Apr 17 '21 at 03:51
  • It's possible I misunderstood your intent. Perhaps I should delete it. – Tim Roberts Apr 17 '21 at 04:00
  • I am expecting `np.asarray(u).T` to return an updated list of `z, theta, dz, dtheta` by the use of the Runge-Kutta method. Meaning, a list of z, list of theta, and so on. – Rafael Gomez Apr 17 '21 at 04:14
  • But `u` is just a list of single values of length t-1. See the last line of the loop in `rk4`? You're appending one value each time through the loop. – Tim Roberts Apr 17 '21 at 04:18
  • Is there a way to get to my goal? I cannot see any form of solution now. – Rafael Gomez Apr 17 '21 at 04:27