0

I am trying to solve this set of differential equations with solve_ivp in Python. Winds are horizontal two dimensional vectors, so I need differential equation solved in y dimension.

enter image description here

This is the error I get:

ValueError: not enough values to unpack (expected 6, got 4)

Is there any other way of solving it? Here is my code:(just projectile motion in x-z plane, distance as function of angle). For the example of two such initial angles (with which we hit approximately the same point) I am trying to calculate the dispersion of hits in the wind according to the point hit in no wind. Thanks for any help or pointers

from random import randrange
from statistics import stdev
from turtle import title
import numpy as np
import scipy as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

    
m = 30
v0 = 301
g = 9.81
ro = 1.28
A = 0.0026
c = 1
k = 0.5 * c * ro * A
B = k / m
V = v0
c_max = 10
c_windx = 0
c_windy = 0

def dSdt(t, S, B, c_windx, c_windy):
    x, vx, y, vy, z, vz = S
    speed = np.hypot(np.hypot((vx + c_windx),(vy + c_windy)), vz)
    return [vx, -B * speed * (vx + c_windx), vy, -B * speed * (vy + c_windy), vz, -B * speed * vz - g]

def distance(angle, B, c_windx, V=v0, t=600):
    v0x = V * np.cos(np.radians(angle))
    v0y = 0
    v0z = V * np.sin(np.radians(angle))
    sol = solve_ivp(dSdt, [0, t], y0 =[0,v0x,0,v0y,0,v0z], t_eval=np.linspace(0,t,10000), args=(B, c_windx,c_windy), atol=1e-7, rtol=1e-4)
    justabove = np.where(np.diff(np.sign(sol.y[2])) < 0)[0][0]
    justunder = justabove + 1
    xloc = (sol.y[0][justabove] + sol.y[0][justunder])/2 
    return xloc

angles = np.linspace(0, 90, 200)
xlocs = np.vectorize(distance)(angles, B=B, c_windx = c_windx)


plt.plot(angles, xlocs)
plt.xlabel('angles []')
plt.ylabel('distance [m]') 
plt.axvline(angles[np.argmax(xlocs)], ls='--', color='r')
plt.text(angles[np.argmax(xlocs)] + 0.2 ,0, round(angles[np.argmax(xlocs)], 2), rotation=90, color='r') 
plt.title('Distance(angle)')
plt.show()

alfa = 15 #its optional

def same_distance_alfa(B,c_windx, V=v0, t=600):
        for a in np.arange(alfa + 5, 90, 0.1):            
            if abs(distance(alfa, B, c_windx, V=v0, t=600) - distance(a, B, c_windx, V=v0, t=600) < 10):
                print((distance(15, B, c_windx, V=v0, t=600) - distance(a, B, c_windx, V=v0, t=600)))
                return a


beta = same_distance_alfa(B)
rand111
  • 15
  • 3
  • Your system is six dimensional, but you only pass four initial conditions to `solve_ivp` with `y0 =[0,v0x,0,v0z]`. – Warren Weckesser Jul 25 '22 at 12:38
  • @WarrenWeckesser thanks. so I fixed my code but there are still some errors and it doesnt work. Is anything else I did wrong? – rand111 Jul 25 '22 at 16:07
  • Your resistance terms in the formulas are wrong. For the velocity `v=v_p+v_w`, the vector resistance is constants times `norm(v)*v`, not just the squares of the separate components. The code has this correctly. – Lutz Lehmann Jul 25 '22 at 18:59
  • If I understand correctly, you want to catch the event(s?) where `y` crosses the zero value in downward direction? The event mechanism of `solve_ivp` can find such points with higher precision than just one step of the bisection method. – Lutz Lehmann Jul 26 '22 at 09:30
  • @LutzLehmann no I want to catch event when `z` crosses the zero value, but also know how did the `x` and `y` change so I know where the projectile lands – rand111 Jul 27 '22 at 10:05
  • The `z` coordinate in your ordering of the components `x,vx,y,vy,z,vz` is `sol.y[4]`. With a defined event, the state at the event roots is recorded in `sol.y_events[0]`. From there you can extract the relevant x and y values. If the event is terminal, the state is also the last in `sol.y`. – Lutz Lehmann Jul 27 '22 at 10:28

0 Answers0