Projectielbaan

[Rik Speybrouck]

Voor wie geinteresseerd is, een verhandeling uit 1977:

https://www.researchgate.net/publicatio ... _the_speed

ik denk dat alle analytische oplossing wel zeer dicht bij elkaar zullen liggen hoor

[Rik Speybrouck]

Wat betreft de horizontale afstand zit er ergens toch een foutje in een excel cel
De fout is ongeveer 2m in het hoogste punt en 4m in het landingspunt
data.png

mijn numerieke benadering komt op 228.51

[wnvl1]

Wat betreft de horizontale afstand zit er ergens toch een foutje in een excel cel
De fout is ongeveer 2m in het hoogste punt en 4m in het landingspunt
data.png

zal aan low angle toepassing liggen denk ik hoe hoger de hoek hoe verder ik kom in vergelijking met numerieke benadering. Jammer dat die formule niet werkt met lambert w functie. Ik heb ook een file met proportionele weerstand waar je op dezelfde manier moeten werk uiteraard met andere formules en daar werkt het perfect. IK zie het in ieder geval niet.

Zoiets oplossen in Python of met een ander wiskundig pakket lijkt mij echt de voorkeur te genieten. Dat is veel gemakkelijker te debuggen. Zo'n excel is veel tijdsintensiever om er fouten uit te halen.

[ukster]

Volgende stap is dan de projectielen zoals een bal op de grond laten botsen en terugstuiten in plaats dat ze nu door de grond gaan. Moet wel kunnen met een drietal extra lijntjes Python code schat ik.

nog een stapje verder....
https://www.maplesoft.com/applications/ ... &view=html
Pythoncode?

[wnvl1]

Ik zal het ook eens proberen in Python, moet even efficiënt kunnen.

[Rik Speybrouck]

Ik zal het ook eens proberen in Python, moet even efficiënt kunnen.

zie jij wat ik fout zou doen bij toepassing van de lambert functie in mijn bericht van eergisteren om de iedeale x te berekenen.

[Xilvo]

Ik zal het ook eens proberen in Python, moet even efficiënt kunnen.

1. voor de x-waarde van het projectiel x_p de y-waarde van de "bodem"-functie y_bodem uitrekenen.
2. Is y_p<y_bodem dan
3. helling bodem bepalen en daaruit normaal op bodem.
4. snelheid projectiel ontbinden in component loodrecht op normaal en component parallel normaal.
5. component parallel aan normaal van teken omkeren.
6. snelheid terugrekenen naar v_x, v_y

[Xilvo]

Dan wordt het zoiets, met de bodemfunctie uit de link van Ukster:

Code: Selecteer alles

def bodem(x):
    return np.sin(x)+(0.2*np.cos(4*x+np.sin(4*x)))-0.2*x+3

def run():
    stuitertel=0
    dt=0.0001
    t=0
    x=0
    y=bodem(x)
    g=9.81
    k=0.0006
    v0=7
    th=40
    thr=np.deg2rad(th)
    vx=v0*np.cos(thr)
    vy=v0*np.sin(thr)    
    print(vx,vy)
    tr=[0]
    xr=[0]
    yr=[y]
    vxr=[]
    vyr=[]
    dbx=0.001
    while stuitertel<5:
        v2=vx**2+vy**2
        Fw=g*k*v2
        Fx=-Fw*vx/np.sqrt(v2)
        Fy=-Fw*vy/np.sqrt(v2)-g
        vx+=dt*Fx
        vy+=dt*Fy
        x+=dt*vx
        y+=dt*vy
        t+=dt
        xr.append(x)
        yr.append(y)        
        tr.append(t)
        vxr.append(vx)
        vyr.append(vy)
        if bodem(x)>y:
            dby=bodem(x+dbx/2)-bodem(x-dbx/2)
            normaal=np.array([-dby,dbx]) # loodrecht op bodem
            normaal/=np.sqrt(np.dot(normaal,normaal)) # eenheidsvector
            v_vect=np.array([vx,vy])
            inprod=np.dot(v_vect,normaal) # inprod. snelheid en normaal
            v_vect-=2*inprod*normaal # component loodrecht bodem omkeren
            vx=v_vect[0]
            vy=v_vect[1]
            stuitertel+=1
                        
    return tr, xr,yr,vxr,vyr

[wnvl1]

Dit is die van mij. Ik heb mij wel eerder gericht op het namaken van de link van ukster in python ipv het originele probleem. De normaal berkenenen doe ik in tegenstelling tot Xilvo door symbolisch af te leiden met sympy ipv numeriek. Als je de code uitvoert zie je een animatie. Ik ben wel vertrokken van gekopieerde code die ik vond voor een vlakke ondergrond, daardoor was het weinig werk. Nu nog in Excel

Code: Selecteer alles

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
import sympy as sym


# Acceleration due to gravity, m.s-2.
g = 9.81
# The maximum x-range of ball's trajectory to plot.
XMAX = 10
YMAX = 6
# The coefficient of restitution for bounces (-v_up/v_down).
cor = 0.99
# The time step for the animation.
dt = 0.005

# Initial position and velocity vectors.
x0, y0 = 2, 4.5
vx0, vy0 = 0, 0

def get_pos(t=0):
    """A generator yielding the ball's position at time t."""
    x, y, vx, vy = x0, y0, vx0, vy0
    while x < XMAX:
        t += dt
        x += vx * dt
        y += vy * dt
        vy -= g * dt
        if y < f(x):
            # bounce!
            y = f(x)
            vector_v = np.array([vx, vy])
            vector_normal = np.array([-f_diff(x), 1])
            vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
            v_reflected = -(1+cor) * np.dot(vector_normal,vector_v) * vector_normal + vector_v
            vx = v_reflected[0]
            vy = v_reflected[1]

        yield x, y

def init():
    """Initialize the animation figure."""
    ax.set_xlim(0, XMAX)
    ax.set_ylim(0, YMAX)
    ax.set_xlabel('$x$ /m')
    ax.set_ylabel('$y$ /m')
    line.set_data(xdata, ydata)
    ball.set_center((x0, y0))
    height_text.set_text(f'Height: {y0:.1f} m')
    return line, ball, height_text

def animate(pos):
    """For each frame, advance the animation to the new position, pos."""
    x, y = pos
    xdata.append(x)
    ydata.append(y)
    line.set_data(xdata, ydata)
    ball.set_center((x, y))
    height_text.set_text(f'Height: {y:.1f} m')
    return line, ball, height_text

# Set up a new Figure, with equal aspect ratio so the ball appears round.
fig, ax = plt.subplots()
ax.set_aspect('equal')

x = sym.Symbol('x')
y_expr = sym.sin(x) + (0.2*sym.cos(4*x+sym.sin(4*x))) - 0.2*x + 3
y_expr_diff = sym.diff(y_expr, x)

f=sym.lambdify(x, y_expr, "numpy")
f_diff=sym.lambdify(x, y_expr_diff, "numpy")

x1 = np.arange(0, 10, 0.01)
y1 = f(x1)

ax.plot(x1, y1, lw=2)

# These are the objects we need to keep track of.
line, = ax.plot([], [], lw=2)
ball = plt.Circle((x0, y0), 0.08)
height_text = ax.text(XMAX*0.5, y0*0.8, f'Height: {y0:.1f} m')
ax.add_patch(ball)
xdata, ydata = [], []

interval = 1000*dt
ani = animation.FuncAnimation(fig, animate, get_pos, blit=True,
                      interval=interval, repeat=False, init_func=init)
plt.show()

[Rik Speybrouck]

Dit is die van mij. Ik heb mij wel eerder gericht op het namaken van de link van ukster in python ipv het originele probleem. De normaal berkenenen doe ik in tegenstelling tot Xilvo door symbolisch af te leiden met sympy ipv numeriek. Als je de code uitvoert zie je een animatie. Ik ben wel vertrokken van gekopieerde code die ik vond voor een vlakke ondergrond, daardoor was het weinig werk. Nu nog in Excel

bal.png

Code: Selecteer alles

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
import sympy as sym


# Acceleration due to gravity, m.s-2.
g = 9.81
# The maximum x-range of ball's trajectory to plot.
XMAX = 10
YMAX = 6
# The coefficient of restitution for bounces (-v_up/v_down).
cor = 0.99
# The time step for the animation.
dt = 0.005

# Initial position and velocity vectors.
x0, y0 = 2, 4.5
vx0, vy0 = 0, 0

def get_pos(t=0):
    """A generator yielding the ball's position at time t."""
    x, y, vx, vy = x0, y0, vx0, vy0
    while x < XMAX:
        t += dt
        x += vx * dt
        y += vy * dt
        vy -= g * dt
        if y < f(x):
            # bounce!
            y = f(x)
            vector_v = np.array([vx, vy])
            vector_normal = np.array([-f_diff(x), 1])
            vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
            v_reflected = -(1+cor) * np.dot(vector_normal,vector_v) * vector_normal + vector_v
            vx = v_reflected[0]
            vy = v_reflected[1]

        yield x, y

def init():
    """Initialize the animation figure."""
    ax.set_xlim(0, XMAX)
    ax.set_ylim(0, YMAX)
    ax.set_xlabel('$x$ /m')
    ax.set_ylabel('$y$ /m')
    line.set_data(xdata, ydata)
    ball.set_center((x0, y0))
    height_text.set_text(f'Height: {y0:.1f} m')
    return line, ball, height_text

def animate(pos):
    """For each frame, advance the animation to the new position, pos."""
    x, y = pos
    xdata.append(x)
    ydata.append(y)
    line.set_data(xdata, ydata)
    ball.set_center((x, y))
    height_text.set_text(f'Height: {y:.1f} m')
    return line, ball, height_text

# Set up a new Figure, with equal aspect ratio so the ball appears round.
fig, ax = plt.subplots()
ax.set_aspect('equal')

x = sym.Symbol('x')
y_expr = sym.sin(x) + (0.2*sym.cos(4*x+sym.sin(4*x))) - 0.2*x + 3
y_expr_diff = sym.diff(y_expr, x)

f=sym.lambdify(x, y_expr, "numpy")
f_diff=sym.lambdify(x, y_expr_diff, "numpy")

x1 = np.arange(0, 10, 0.01)
y1 = f(x1)

ax.plot(x1, y1, lw=2)

# These are the objects we need to keep track of.
line, = ax.plot([], [], lw=2)
ball = plt.Circle((x0, y0), 0.08)
height_text = ax.text(XMAX*0.5, y0*0.8, f'Height: {y0:.1f} m')
ax.add_patch(ball)
xdata, ydata = [], []

interval = 1000*dt
ani = animation.FuncAnimation(fig, animate, get_pos, blit=True,
                      interval=interval, repeat=False, init_func=init)
plt.show()

De excel file is wel gemaakt op basis van een zelf uitgewerkte differentiaal vergelijking type low angle gewoon met pen en papier. Behoudens een kleine aanpassing inzake beginhoogte die ik heb doorgegeven aan ukster is hij perfect. En zelfs onder een hoek van 45° is de afwijking in % nog niet zo heel hoog.( uitwerking werd verleden week on line gezet)

[OOOVincentOOO]

@wnvl1,

Nice , erg leerzaam met animation en sympy. Vaak worden alleen plaatjes getoond zonder de code/methode of half (voor mij niet interessant, ik wil onder de motorkap kijken!). In dit topic soms dus wel mogelijk te kijken hoe iedereen zijn oplossing vind.

Echter er klopt iets niet met een offset wanneer ik bij mij code lopen laat (zie plaatje). Heb mij verder niet verdiept in topic wat de oorzaak zou kunnen zijn. Heb kort rond gespeelt maar niets gevonden.

nb.
In de code header voeg ik het volgende toe om een pop scherm te krijgen qt5. Lijkt mij onwaarschijnlijk de reden te zijn. Maar men weet nooit.

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

[wnvl1]

@Rik Speybrouck

Ik wilde overigens niet bedoelen dat dat niet lukt in Excel. Ook de berekeningen met deze ondergrond krijg je nog wel zonder al te veel problemen gedaan binnen Excel. Ik wilde in een eerder berichten aangeven dat het moeilijker is om te debuggen omdat het minder overzichtelijk is. Als ik de code van Xilvo zie dan heb ik in een oogopslag door hoe ze werkt hoewel ik ze niet zelf geschreven hen, hetzelfde bij de Maple code van Ukster, dat gaat bij een Excel niet lukken.

[wnvl1]

nb.
In de code header voeg ik het volgende toe om een pop scherm te krijgen qt5. Lijkt mij onwaarschijnlijk de reden te zijn. Maar men weet nooit.

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

Ik gebruik Python 3.8 in Pycharm. Bij mij opent dat vanzelf een apart window waarin de achtergrond en de animatie verschijnt. Ik kan het probleem niet reproduceren. Is wel een vervelende fout.

[wnvl1]

Ter aanvulling als ik de configure subplots (derde van recht bovenaan op screenshot OOOVincentOOO) gebruik tijdens de animatie dan heb ik het probleem ook. Het window gewoon resizen geeft geen probleem.

[OOOVincentOOO]

@wnvl1,

Wil verder niet teveel off topic. Ik gebruik Spyder Python 3.8.

Als ik XMAX en YMAX aanpas tot grafiek bereik welke 10 bij 5 is krijg ik iets betere resulaten maar nog niet helemaal.

Code: Selecteer alles

XMAX =10
YMAX = 5

Verder forceer je de assen *equal* en begrjip ik ook:

Code: Selecteer alles

ax.set_aspect('equal')

Dat betekend dat een van de assen een compromis moet sluiten en zich aanpassen in limieten. Heb geprobeerd met:

Code: Selecteer alles

YMAX=ax.get_ylim()

Word iets beter maar nog niet helemaal.

Ik laat het hierbij betreffende praten over Python in dit topic! Wil niet topic verstoren. Heeft weinig te maken met Projectiel baan nu.

EDIT:
Vergeet commentaar hierboven!!
De code werkt goed. Ik had het pop-up grafiek venster erg klein (default afmetingen). Dan gaat er iets mis met schalering!!!
Misschien hebt een fig_size aan tegeven bij init ax.