0

the equation is:

d^2 r/dt^2 = -c/m (dr/dt)+g

where r is the position of the projectile, c is the drag coefficient, m is the mass of the projectile and g is the acceleration due to gravity.

Assuming only two-dimensions in component form this, of course, reads,

d^2 X/dt^2 = -c(dX/dt)= U
d^2 Y/dt^2 = -c/m(dY/dt)+g

if we employ the methodology described above and define explicitly the velocities in X AND Y as,

U = dX/dt

and

V = dX/dt

then the entire coupled system of equations is,

dU/dt= -c/m(U)
dV/dt= - c/m(V)+g
dX/dt= U
dY/dt = V

The parameters for this system of ODEs are c = 0.5 kgs^−1, m = 2kg and g = −9.81 ms^−2.

initialising the variables as (U0, V0, X0, Y0) = (173, 100, 0, 0) which launches the projectile from the origin at an angle of ∼ 30 degrees from the horizontal.

how to I write a new function in python using rk4 (I want to know how to code this) that implements the system of four ODEs above that solve the 2D projectile motion problem....? please help I am very new to ODEs AND CODING. THANKS

I have got the following so far...and its not working and I really don't know what to do for this specific problem, I am meant to obtain a projectile graph as well...can someone please improve my code pls thanks

import numpy as np
import matplotlib.pyplot as plt


def projectileMotion_V(t, M, g, c):
   return -c/M * V0 + g 


def projectileMotion_U(t, c, M):
   return -c/M * U0

V0 = 100         
U0 = 173
ang = 30.0      
c = 0.5       
dt = 0.1  
M = 2.0         
g = -9.81 
h = 0.1

t = [0]                         
x = [0]                         
y = [0]
vx = [V0*np.cos((ang*np.pi)/180)]  
vy = [U0*np.sin((ang*np.pi)/180)]
ax = [-(c*V0*np.cos((ang*np.pi)/180))/M]          
ay = [g-(c*U0*np.sin((ang*np.pi)/180))/M]



def solveODEsWithR4Method(t, x, y, vx, vy, ax, ay):
   t.append(t[0]+dt)                
   vx.append(vx[0]+dt*ax[0])  
   vy.append(vy[0]+dt*ay[0])
   x.append(x[0]+dt*vx[0])    
   y.append(y[0]+dt*vy[0])    
   vel = np.sqrt(vx[0+1]**2 + vy[0+1]**2)   
   drag = c*vel                                    
   ax.append(-(drag*np.cos(ang/180*np.pi))/M)     
   ay.append(-g-(drag*np.sin(ang/180*np.pi)/M)) 
return -c/M * V0 + g 


fig,ax = plt.subplots()
ax.plot(t, M, g, c)
plt.show()
E_net4
  • 27,810
  • 13
  • 101
  • 139
TIVONA
  • 19
  • 1
  • 1
    Hi there, and welcome- What has your 'net search turned up? For example, I see: https://www.codeproject.com/tips/792927/fourth-order-runge-kutta-method-in-python , and https://www.youtube.com/watch?v=IOkwWYaZbck There are two things here, coding, and solve the ODEs. You might want to start with coding something simple, like the solution to the quadratic equation, and then a simple 1D ODE, to get your 'python legs', and then tackle this problem. – asylumax Apr 20 '20 at 23:29
  • Air resistance is quadratic in the speed, you are using water resistance (but with the wrong coefficient), see https://stackoverflow.com/a/35009733/3088138 for how that computation should look like. – Lutz Lehmann Apr 21 '20 at 06:05
  • This is an inhomogeneous system of linear differential equations with constant coefficients. You can solve it completely and obtain the exact solutions, so there is no need for numerical integration. – Futurologist Apr 24 '20 at 14:07

2 Answers2

3

As nothing else happens in these times of a pandemic of delusions about phantasmic killer viruses...

Please do not call it air resistance, that is specifically k*|v|*v. As for this to be a force the coefficient k would need to have units kg/m, which is not what you are given, your formulas for the resistance are probably correct. Call it "medium resistance" instead, water resistance would behave that way.

Then code the acceleration

c = 0.5; m = 2; g = -9.81;
def motion(x,v):
    x,y,vx,vy = v
    return np.array([vx,vy, -c/m * vx + g, -c/m * vy ])

Copy the RK4 code from somewhere that is geared towards vector states 

def RK4step(f,u,dt):
    k1 = dt*f(u)
    k2 = dt*f(u+0.5*k1)
    k3 = dt*f(u+0.5*k2)
    k4 = dt*f(u+k3)
    return u + (k1+2*k2+2*k3+k4)/6

def RK4integrate(f, u0, tspan):
    u = np.zeros([len(tspan),len(u0)])
    u[0,:]=u0
    for k in range(1, len(tspan)):
        u[k,:] = RK4step(f, u[k-1], tspan[k]-tspan[k-1])
    return u

and apply both codes together to compute a trajectory

dt = .1
t = np.arange(0,10,dt)
u0 = np.array([0, 0, 173, 100])

sol_RK4 = RK4integrate(motion, u0, t)
x,y,vx,vy = sol_RK4.T
plt.plot(x,y)
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
0

Set v = dr/dt. Then the equations become:

dv/dt = - (c/m) * v  +  g

Multiply both sides of the equation with exp((c/m)*t):

exp((c/m)*t) * dv/dt = - exp((c/m)*t) * (c/m) * v  +  exp((c/m)*t) * g

Move the first column of terms from the right-hand side to the left one:

exp((c/m)*t) * dv/dt + exp((c/m)*t) * (c/m) * v  =  exp((c/m)*t) * g

Then by the product rule of differentiation, applied to exp((c/m)*t) * v, one gets 

d/dt( exp((c/m)*t) * v )  =  exp((c/m)*t) * g

d/dt( exp((c/m)*t) * v )  =  d/dt( (m/c) * exp((c/m)*t) * g )

Now you can integrate both sides with respect to t and there exists a vector u0 such that

exp((c/m)*t) * v  =  u0  +  (m/c) * exp((c/m)*t) * g

Multiply both sides by exp(-(c/m)*t):

v  =  exp(-(c/m)*t) * u0  +  (m/c) * g

If the initial velocity is v0, then you have to set u0 = v0 - (m/c)*g so the final formula for the velocity is:

v  =  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * g

Integrate v one more time with respect to t and you get the formula for the position:

r  = r0 - (m/c) * ( v0 + (m/c)*g )  - (m/c) * exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * t * g

where r0 is the initial position.

So the final formulas for the solution are

r  = r0 - v0 + (m/c)*g  +  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * t * g

v  =  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * g
Futurologist
  • 1,874
  • 2
  • 7
  • 9